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Preface 



This volume contains the proceedings of the Fourth Workshop on Flyhrid Sy- 
stems: Computation and Control (HSCC 2001) held in Rome, Italy on March 
28-30, 2001. The Workshop on Hybrid Systems attracts researchers from indu- 
stry and academia interested in modeling, analysis, synthesis, and implementa- 
tion of dynamic and reactive systems involving both discrete (integer, logical, 
symbolic) and continuous behaviors. It is a forum for the discussion of the la- 
test developments in all aspects of hybrid systems, including formal models and 
computational representations, algorithms and heuristics, computational tools, 
and new challenging applications. 

The Fourth HSCC International Workshop continues the series of workshops 
held in Grenoble, France (HART’97), Berkeley, California, USA (HSCC’98), Nij- 
megen, The Netherlands (HSCC’99), and Pittsburgh, Pennsylvania, USA (HSCC 
2000). Proceedings of these workshops have been published in the Lecture Notes 
in Computer Science (LNCS) series by Springer- Verlag. 

In line with the beautiful work that led to the design of the palace in which 
the workshop was held, Palazzo Lancellotti in Rome, resulting from the colla- 
boration of many artists and architects of different backgrounds, the challenge 
faced by the hybrid system community is to harmonize and extract the best from 
two main research areas: computer science and control theory. Terminology, ma- 
thematical tools, and abstractions are different, problems considered relevant by 
one community may be considered trivial by the other, yet it is this very diffe- 
rence that may bring new vistas to traditional research fields to escape the trap 
of routine. The steering committee of the workshop series has been appointed to 
guide the directions of the research in troubled water balancing the membership 
among computer scientists, control theorists, and application experts. The tech- 
nical program committee has been assembled following the same principle. The 
committee has done a wonderful job in reviewing and discussing 82 submissions 
(a record number since the inception of the workshop series). All requested re- 
views were received (a world-wide record among all workshops!). After extended 
and, at times, intense discussions, 36 papers were selected for presentation at the 
workshop and publication in this volume. While the technical quality of the pa- 
pers is excellent, we cannot underestimate the preponderance of control theory 
papers and the scarcity of application papers. The theory papers are mainly di- 
rected at the consolidation of the foundations of the field, a hardly unexpected 
outcome in an area that is approaching a new level of maturity. However, the 
lack of relevant application papers is somewhat worrisome. For this reason, we 
preferred to give emphasis to applications in the invited papers to the workshop: 
Manfred Morari (ETH Zurich), Costas Pantelides (Imperial College), and Janos 
Sztipanovits (Vanderbilt University) are all well known for their work in hybrid 
system applications and in embedded-system design. In addition, we included 
in the workshop a panel on applications of hybrid systems. The participants to 
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the panel addressed the challenges of using a richly expressive theory, being, as 
such, relatively poor in computationally affordable synthesis and analysis tools, 
to yield relevant results in the real-life world. They also addressed the issue of 
merging knowledge about tools and methods in control and computer science 
so that we may avoid the risk of re-inventing in one field results that are well 
known in the other. 

We believe that embedded systems will be the main application vehicle for 
our technology and as such deserve particular attention. Embedded systems will 
also be the main application domain for electronics in general. Since embedded 
systems require design methods that guarantee correct and efficient behavior in 
harsh environments, a strong theoretical approach to synthesis and verification 
is badly needed. They are hybrid in nature: continuous and discrete mix freely 
in a variety of application domains. Software and control will play a dominant 
role. Hence, we believe that our community will be an important constituency 
in founding the field of embedded system theory and design. 

We wish to thank the organizations (PARADES, Progetto Finalizzato Ma- 
dess II, Consiglio Nazionale delle Ricerche, Army Research Office, National 
Science Foundation) that financially supported the workshop. Moreover, we ack- 
nowledge the contribution of Magneti-Marelli, an automotive electronics com- 
pany that has put to good use hybrid system technology in its products. In 
particular, the support and continuous encouragement of Dr. Daniele Pecchini, 
President and General Manager of Magneti Marelli Powertrain Division, is ack- 
nowledged. We thank Prof. Richard Gerber for letting us use START, his soft- 
ware conference manager. 

The final remark is dedicated to the Organizing Gommittee, whose members 
spent long hours making sure everything was correctly handled, from call for pa- 
pers to hotel information, and paper submission. In particular, Andrea Balluchi 
and Luca Benvenuti have spent an inordinate amount of time coping with the 
software, trying to keep all the web material in synch and making sure authors 
submitted the correct versions of their papers and the appropriate documents 
that E-conomy bureaucracy imposes on us. 
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Based on work jointly with 

Alberto Bemporad, Francesco Borrelli, Francesco Cuzzola, Tobias Geyer, 
Domenico Mignone, Fabio Torrisi and Giancarlo Ferrari Trecate 



Abstract. We envision the role of control to expand rapidly in two di- 
rections. It will impact novel application areas, which have yet to benefit 
from the power of feedback, and, as an embedded technology, control 
will extend its reach far beyond the traditional narrow concept to in- 
clude higher level functions of operation. Our research program is built 
on this vision. Eventually, these ideas should also radically change what 
is taught in our class rooms, so that our students can transfer these 
techniques to industry effectively and reap its benefits. 

In all control applications the actual control algorithm is just one tiny 
part of the overall system designed to ensure safe, reliable and economi- 
cal operation. Success or failure of "operation” are attributable at least 
as much to "the rest” as to the control algorithm itself. At the lowest 
level the control algorithm is endowed with functionality to deal with op- 
erating constraints and to switch smoothly between different operating 
regimes. At the highest levels the control algorithm may be embedded 
in a scheduling system or even an Enterprise Resource Planning (ERP) 
system. At all levels this embedding creates a heterogeneous system com- 
prised of many interacting subsystems, typically referred to as a hybrid 
system. 

The integration should eventually lead to a safer, smoother, more re- 
sponsive and more competitive functioning of the entire system or or- 
ganization. About three years ago we embarked on a major research 
program toward this goal. Its objective is the development of new theo- 
retical tools to model, analyze, simulate and control such large complex 
hybrid systems involving continuous and discrete states, whose behavior 
is governed by dynamics, logical statements and constraints. In this talk 
we will summarize the highlights and try to put them in perspective. 
Modeling and Simulation: The models should facilitate the analysis and, 
at the same time, capture the complex behavior, that hybrid systems are 
known to exhibit. Based on these considerations, we introduced a discrete 
time description, combining linear dynamics with Boolean variables. This 
mixed logical dynamical form (MLD) form is capable to model a broad 
class of systems arising in many applications from the automotive, air- 
craft, chemical and information technology helds. Supply chains used in 
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business models can be conveniently modeled as MLD systems as well. 
We defined a new modelling language (HYSDEL) and wrote a compiler 
to assist the user in the formulation of MLD models. 

Controller Synthesis: For controller synthesis we formulate a finite hori- 
zon optimal control problem and apply the result in a moving horizon 
fashion. For MLD models the optimization problem is a mixed-integer 
linear program (MILP) which must be solved in real time at each sam- 
pling time. We have proven that the resulting state feedback control law 
is piece-wise linear over a polyhedral partition of the state space. As an 
alternative to on-line optimization, we can determine this control law 
explicitly by solving a multi-parametric MILP. 

State Estimation and Fault Detection: For application of the described 
control law the system states must be known. Estimation of the states 
of an MLD system is a complex nonlinear filtering problem. We have 
defined a moving horizon estimator, where at each time step a mixed in- 
teger quadratic program must be solved to arrive at the state estimates. 
We have proven the convergence of the estimator if certain observability 
properties are satisfied. Complex fault situations can be modeled accu- 
rately in the MLD framework. Fault detection is another application of 
the new estimator. 
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Abstract. Most processes of practical interest are hybrid in natnre, ex- 
hibiting both continuous and discrete characteristics. In many cases, the 
hybrid behaviour is a result of intrinsic physical phenomena that lead 
to (practically) instantaneous events such as the appearance and disap- 
pearance of thermodynamic phases, changes in flow regimes, equipment 
failures etc. All such events effect qualitative changes in the underlying 
continuous dynamics, thereby leading to hybrid macroscopic behaviour. 
In other cases, the hybrid nature arises from external discrete actions 
imposed on the process by its control system. For example, the latter 
may apply quantisation to convert continuous process measurements into 
discrete ones and/or continuous control outputs into discrete actions. 
Hybrid processes and hybrid controllers, and their combination, can be 
modelled in terms of State- Transition Networks (STNs). The system be- 
haviour in each state is described by a different set of continuous equa- 
tions (typically a mixed system of partial and/or ordinary differential 
and algebraic equations). At any particular time during its operation, 
the system is in exactly one such state. An instantaneous transition to a 
different state may take place if a certain logical condition becomes true. 
Each transition is also characterised by a set of continuous relations that 
determine unique values for the system variables immediately following 
the transition in terms of their values immediately preceding it. 

In this presentation, we consider mathematical formulations and tech- 
niques for the optimisation of hybrid systems described by STNs. This 
generally seeks to determine the time variation of a set of controls and/or 
the values of a set of time-invariant parameters that optimise some as- 
pect of the dynamic behaviour of the system. The time horizon of interest 
may be fixed or variable, subject to specified lower and upper bounds. 
The equations that determine the system behaviour in each state may 
be augmented with additional inequality constraints imposing certain 
restrictions (related to safety or operability) on the acceptable system 
trajectories. The objective function to be minimised or maximised is usu- 
ally a combination of fixed contributions (depending on the values of the 
time-invariant parameters) and variable contributions (depending on the 
system trajectory, including the variation of the controls). 

As an illustration, we start with simple linear systems operating in the 
discrete time domain, possibly involving uncertain parameters. We then 
proceed to consider the more complex problem of the optimisation of 
nonlinear hybrid systems operating in the continuous time domain. 
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Embedded Software and Systems 
Challenges and Approaches 



Janos Sztipanovits 
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Vanderbilt University, Nashville, TN, USA 



Abstract. One of the most pervasive applications of computing is infor- 
mation processing tightly integrated with physical processes. Embedded 
computing rapidly takes over the role of being a universal integrator for 
physical systems. This trend is based on a fundamental technical reason: 
digital information processing is uniquely suitable for controlling and 
implementing complex interactions among physical system components. 
The expanding integration role of computing challenges the state-of-the- 
art in both system and software design. Eirst, the traditional separation 
of related design disciplines is not maintainable. Predictability of the 
design requires integrated modeling and analysis of physical processes 
and information processing. Second, the narrow focus of current software 
technology on functional composition is not sufficient. Essential physi- 
cal properties of embedded computing systems, such as timing, noise or 
fault behavior, cut across functional boundaries, which makes software 
design and implementation extremely hard and expensive. Third, design 
technologies, which are based on the modeling and analysis of systems 
with static structure, are becoming inadequate. Although networked em- 
bedded computing combined with inexpensive MEMS-based sensors and 
actuators make the construction of large physical systems with continu- 
ously changing structure and physical interactions feasible, their design 
is an open challenge. 

The first part of the talk provides an overview of the unique challenges 
and new research directions in embedded system and software design. 
The second part of the talk describes the Model-Integrated Computing 
(MIC) approach to address some of these challenges. Using the design 
of structurally adaptive embedded processing systems as example, the 
following three topics will be covered: 

1. Methods and tools for the specification and construction of multiple- 
view, domain-specific modeling languages and integrated design en- 
vironments. The MIC approach is based on the application of meta- 
modeling, meta-programmable modeling tools and model translators 
that form the foundation for composable design environments. 

2. Automated synthesis of processing architectures satisfying multiple 
functional and physical constraints. The method described is based 
on symbolic constraint satisfaction. 

3. Application of generative programming techniques with special em- 
phasis on model-based software generators. 
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Hybrid Systems Applications: An Oxymoron? 



Moderator: 

Alberto Sangiovanni-Vincentelli^’^ 

Partecipants: 

Thomas A. Henzinger^, Bruce H. Krogh^, Oded Maler"^, 

Manfred Morari®, Costas C. Pantelides®’”^, George J. Pappas®, 

Tunc Simsec^, Janos Sztipanovits®, Stavros Tripakis®’^® 

^ EECS Department, University of California at Berkeley, CA, USA 
^ PARADES, Roma, Italy. 

® ECE Department, Carnegie Mellon University, Pittsburgh, PA, USA 
VERIMAG, Centre Equation, Gieres, France 
® Automatic Control Laboratory, ETH, Zurich, Switzerland 
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^ Process Systems Enterprise Ltd., London, United Kingdom 
® EE Department, University of Pennsylvania, PA, USA 
® EECS Department, Vanderbilt University, Nashville, TN, USA 
CALIFORNIA PATH, University of California at Berkeley, CA, USA 

Hybrid systems are richly expressive models for a large variety of potential ap- 
plications. However, being so rich as to include continuous nonlinear dynamical 
systems, discrete-event systems and other models of computation (finite-state 
machines and data flow come to mind here), they are not amenable to com- 
putationally attractive techniques for synthesis and analysis and present hard 
numerical problems to simulation. Hence, applying the methods typical of this 
technology requires non trivial amount of approximation and abstraction. And 
approximation and abstraction are effective only if the domain of application 
is deeply understood. Thus, significant applications of hybrid systems require a 
great deal of work both to select the right abstraction level and to derive algo- 
rithms that exploit the particularities of the domain of application. In addition, 
one needs to motivate and document convincingly why using hybrid systems can 
yield better results than other techniques. In this respect, there has been an on- 
going debate as to what constitutes a meaningful result in applications: on one 
hand, novel languages for describing hybrid systems and capturing their prop- 
erties may be considered sophomoric exercises by experts in languages, on the 
other, formal verification tools that in general can handle small systems may be 
seen as toys for who is trying to tame entire chemical plants. On the simulation 
front, how to deal with discontinuities of trajectories is a major issue. Numerical 
analysts have been looking at these problems only recently and with a great deal 
of skepticism as to what can be proven rigorously. Hybrid system researchers are 
now getting seriously in the simulation arena exploiting what has been done in 
the numerical analysis arena. 

The goal of the panel is to bring experts from the two reference communi- 
ties of hybrid systems (computer science and control) to debate whether hybrid 
system applications can indeed be compelling and what can be done to prevent 
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naive work on both sides when straddling across competence domains. Simula- 
tion and verification in general will also be discussed in the frame of the work 
done in numerical analysis. Predicting the outcome of the panel, we would like 
to end with a positive note: hybrid system technology is relevant to important 
applications but it has to be handled with great care and pushing the cart all in 
the same directions will give the hybrid system community the relevance it has 
the right of aspiring to. 
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Abstract. An approach to estimation for a class of hybrid discrete-time 
linear systems using Luenberger observers is presented. The proposed 
Luenberger observer for such a kind of systems relies on the switching 
among different gains. Convergence conditions have been found to ensure 
the stability of the error dynamics and the related gains may be selected 
by solving a set of linear matrix inequalities (LMIs). Moreover, this ob- 
server may be improved by suitably updating the estimate nsing the last 
measures. This update enables one to reduce the norm of the estimation 
error and is based on the so-called projection method. Simulation results 
are reported to show the effectiveness of these methods in the estimation 
for hybrid discrete-time linear systems. 



1 Introduction 

Hybrid systems have recently gained a great attention and the research in this 
area has been devoted more to control problems. In this work, the subject is the 
state estimation for a class of hybrid systems described by switching discrete-time 
linear equations. Switching systems are well-suited to dealing with applications 
like, for example, gain scheduling, reconfigurable control, and fault diagnosis, 
which enable one to point out the importance of constructing observers for such 
systems. 

The problem of estimating the state of a switching system was originally 
stated in fP . Later on, a lot of researches investigated the issues related to such 
a problem in a probabilistic framework, i.e., supposing that the transitions occur 
according to a model described by a first-order Markov chain. Difficulties may 
arise in the solution of optimal Bayesian estimation problems and the interested 
reader is referred, among others, to |2] and |3|. Another relevant topic concerns 
the so-called multi-model estimation. Such a subject is quite vast and involves 
many application-oriented problems (for an introduction, see 0). Summing-up, 
all the above-mentioned approaches rely on a stochastic setting and the switching 
event is supposed unknown. Here, we focus on the problem of estimating the 
state of a switching system by assuming to know both time and mode of the 
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switching. Nevertheless, also in this context, the problem remains hard to solve 
for the difficulties of both guaranteeing the stability of the estimation error and 
devising a suitable, efficient observer design procedure. It is worth noting that 
we will make no assumption on the probabilistic description of the system mode 
transitions to derive the stability results of the proposed estimation methods. 

Gain switching observers for continuous-time nonlinear systems have been 
considered in where stable switching laws are searched for with different 
Lyapunov functions for each gain. A different approach based on coprime factor- 
ization is proposed in |E| to construct an observer for switching continuous-time 
linear systems. In the present paper, the goal is to find an estimator with a 
stable estimation error in the presence of any switching in a given finite set of 
admissible system modes. Such a problem turns out to be more difficult than 
the standard design of Luenberger observers for time-invariant linear systems. 
In this case, a Luenberger observer provides a convergent error dynamics if and 
only if the gain is chosen such that the poles of the error dynamics are in the 
strictly stable region. This condition is not sufficient to ensure the stability of 
the estimation error for a switching linear system. The gain selection of a switch- 
ing observer is nontrivial as it involves the typical stability issues of the hybrid 
systems (for an introduction, see 0, [5|, and 0). In our case, the solution of 
this problem has been addressed by seeking a common Lyapunov function. This, 
in turn, can be reduced to the fulfillment of linear matrix inequalities (LMIs), 
which allow one to easily obtain a solution in a computationally feasible way. 

An improvement to this Luenberger observer has been made by applying a 
projection method imiTi to update the current estimate using the last mea- 
sures. The resulting estimator exhibits a stable error dynamics if the same LMI 
relationships found for the first estimator are satisfied. In addition, the new ob- 
server results in higher performance, as this update provides a reduction of the 
estimation error. 

The paper is organized as follows. Section |21 is devoted to the problem of 
constructing a Luenberger observer for switching discrete-time linear systems, 
with a particular emphasis on the stability of the error dynamics and on the 
development of an LMI approach to synthesize such observers. In Section 01 
a modified Luenberger observer with the related stability analysis is proposed 
that enables one to estimate the state of the system using also the last available 
measures. Simulation results are illustrated in Section 0to show the performance 
of the proposed estimation methods. The conclusions are drawn in Section 0 

2 Switching Observers for Discrete-Time Linear Systems 

Consider the discrete-time linear system 

r Xt+i = A^(t) Xt + t = 0, 1, . . . (1) 

where G R” is the state vector, Uf. G is the input vector, G R™ 
is the measure vector, and a : N — > {l,2,...,fc} is a function that maps 
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the index time stage into an index set {1, 2, . . . , fc} . Each of the indices cor- 
responds to a different model of the system and measurement equations, i.e., 

^cr(t) € A = {Ai, A 2 , . . ■ , Ak} , £ B = {Bi, B 2 , ■ ■ ■ , Bk} , C 

C = {Ci,C 2 ,...,Ck}, where Ai G B, G R'^^p , and C, G for 

i = 1,2, . . . ,k ■ We assume that the matrices Ci G i = 1, 2, . . . , fc , are 

of full rank m < n and the output of the function cr(-) is known at time t. 
Anyway, it is worth noting that the matrices B^(^t-^ and C„(t) with time-varying 
dimensions in the number of columns and rows (i.e., of p and m, respectively) 
are allowed. A switching observer for (0 is the following: 

At “k (y_t ^a(t) At) ) ^ 0,1,... (2) 

where Iq = ig is chosen “a priori” and is the observer gain at the time 

t , Ba{t) € £ = {Lx,L 2 , . . . , Lk] , and Li G , i = 1, 2, . . . , fc . A pictorial 

representation of such an observer is shown in Fig. ^ 




Fig. 1. Scheme of a switching observer. 



Note that these gains may change in such a way that the dimension m will 
vary over time due, for example, to a variable number of available measures at 
time t. The dynamics of the estimation error (i.e., ^t= At ~ At) behaves like a 
switching dynamic system, thus a common Lyapunov function is searched for to 
ensure stability. Now, we can state the following theorem. 



Theorem 1. Consider the system (Cj) and assume that the pairs (Ai,Ci) , i = 
1,2, ... ,k, are observable. If there exists a symmetric positive definite matrix P 
as the solution of the algebraic Lyapunov inequalities 

{A,-L,Ci)'^ P{A,-L,Ci)-P <D , i = l,2,...,k (3) 

then the observer involves an estimation error asymptotically convergent to 



zero. 
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Proof of Theorem ^ Let us prove the result stated in Theorem D The 
error dynamics may be computed by means of equations (P) and m 



If we consider the Lyapunov function Vt = ef P , where P is a symmetric 
positive definite matrix, we obtain Vt+i < Vt , Ve^ G R" , if 



and then (0| may be easily derived. 

■ 

It is important to recall that the assumption on the observabilty of the pairs 
(Ai, Ci) , J = 1, 2, . . . , fc, is necessary to guarantee that each inequality in @ may 
admit a solution for a given positive definite matrix P , but the existence of a 
common P satisfying all the inequalities is required to ensure stability. As it is 
difficult to find a common Lyapunov function once the gains Li , i = 1,2, . . . ,k, 
have been selected, we will try to find the gains and the positive definite matrix 
P simultaneously. Thus, the goal is to solve the following problem. 

Problem 1. Find Li , i = 1,2, . . . ,k, such that there exists a symmetric positive 
definite matrix P solving the Lyapunov inequalities 



The above problem may be reduced to a simpler form that is well-suited 
to being solved by an LMI method. To this end, let us consider the following 
lemma. 

Lemma 1. Given a symmetrie positive definite matrix P , an inequality 



(,Afj(^t) ^a{t) ^cr(t)) §Lt : I 0, 1, . . . 



(^o-(t) ~ ^cr(t) P (^o-(t) ~ ^cr(t) C'o-(t)) ~ P < 0 , t — 0,1, . . . 



{A,-L,af P{A,-L,a)-P <0 , i = l,2,...,k . (4) 



□ 



{A, - L, P {A, -L,a)-P<0 



( 5 ) 



is equivalent to 




( 6 ) 



where Li = P ^ fy , i = 1,2, . . . ,k. 



□ 



Proof of Lemma ID Let us recall the well-known Schur complement 

Q-SR-^S'^ >0 o Q > 0, R- Q-'^ S > 0 
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where i? , S , and Q are matrices of appropriate dimensions. If we apply this 
result by taking Q = P , S = P Ai —YiCi , and R = P , we can easily verify 
that 0 gives © if Li = y, . 



To sum up, the solution of Problem ^can be obtained by solving the following 
LMI problem. 



Problem 2. Find P > 0 and Yi, i = 1,2,. 

/ P PA,-Y,CA 

{(PA-Y^af P J 



■ ■ ,k, such that 



>0 , i=l,2,.. 



.,k 



and take the observer gains Li = P ^ Yi . 



(7) 



□ 

ProblemQis simpler than Problem Q as the former is linear in the unknown 
parameters, whereas the latter is quadratic at the first glance. Moreover, the 
formulation of Problem E| fits the so-called LMI framework which enables 
one to solve it by means of convex programming algorithms. Efficient numerical 
methods for convex optimization are available, and the reader is referred to H2j 
for an introduction on this subject. 



3 An Enhanced Projection-Based Luenberger Observer 

The Luenberger observer (0 provides an estimate of the state at time t -I- 1 
using the measures available at time t by means of . As a matter of fact, we 

aim at determining the estimate using also like a standard Kalman 

filter. To this end, a method is proposed and consists in updating the estimate 
given by the Luenberger observer © by means of the projection method | n 1 1 I j . 
which allows one to take into account the last measures. More specifically, this 
estimation method is performed as follows: 

J Lt+i = Lt + 7?o-(t) Mt + {y_t ~ 1 t = 0 , 1 , . . . 

|Lt+i=Lt+i+L (c'(r(i-i-i)L {Vt+i ~ Lt+i) 

( 8 ) 

where ^ is chosen “a priori”, Pcr(t) is the observer gain at time t, 

i.e., Pcr(t) & L = {Li, L 2 , ■ ■ ■ , Lk} , and P is a positive definite matrix. The 
update that enables one to derive the new estimate using and the 

last measure vector is based on a simple geometrical idea we will illustrate 
in the following. 

For the sake of notational simplicity, let assume x^_^_i is measured at time 
t-l-1 by means of = C x^_^_^ (i.e., C is used instead of ). Moreover, 
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we regard ij+i as an “a priori” estimate of at time t + 1 and want to 
determine a new estimate such that ||e^]^|| < ||ej_|_j^||, where = 

Xt+i ~^+i ■ The state space can be decomposed into two orthogonal subspaces, 
like, for example, the null space of C (i.e., N {C) = { x G R” : C a: = 0 } ) and its 

orthogonal space N using the scalar product < x,z >p= ^ P z, x,zG 

K" (this scalar product is well-defined as the matrix P is positive definite). 
If P is taken equal to the identity matrix, it is easy to verify that N (C)'*" is 
R (C^) (i.e., the space spanned by linear combinations of the columns of the 
matrix C '^ ) . 




Fig. 2. Sketch to explain the projection method {C replaces C'o-(t+i)). 



The decomposition can be accomplished by means of the subspaces given by 
R{P~^ C^) and its orthogonal complement, instead of N (C')"'" and N (C) . 
The reason for using this subspace decomposition concerns the stability of 
the estimation error as it will be clarified in the following. Fig. 0 pro- 
vides a meaningful geometrical interpretation of the projection method and 
enables one to illustrate the rationale for the proposed approach. As can 
be noticed in Fig. El the projection of on R(P~^ C'^) is equal to 

P-^C^ (CP-^C^)~^ i.e., P-^C^ (CP-^C'^)~^ . Note that the 

projection matrix P~^C^ (CP~^C'^) ^ is well-defined as the matrices Ci € 
z = 1,2, ...,fc (i.e., C in this case) have been assumed of full rank 
m < n . In practice, the estimate of x^_^^ is obtained by projecting x^J^^ on the 
subspace corresponding to the new measure , which provides a new esti- 
mate x^pi such that the corresponding estimation error is smaller than 

that of the previous error, i.e., ||e)'l|_j^|| < ||e(_|_]^|| . A pictorial representation of 
the observer Q is shown in Fig. 0 

As far as it concerns the stability of estimation error associated to (0, we 
can state the following result. 
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Fig. 3. Scheme of a switching projection-based observer. 



Theorem 2. Consider the system ^ and assume that the pairs (Ai,Ci) , i = 
1,2, ... ,k, are observable. If there exists a symmetric positive definite matrix P 
as the solution of the algebraic Lyapunov inequalities 



{A,-L,C^f P{A,-L,Ci)-P <Q , i = l,2,...,k (9) 

then the estimator 0 involves an estimation error asymptotically convergent to 
zero. 



□ 

Proof of Theorem 121 The estimation error before the projection update is 
given by 



Ht+l — (^o-(t) C'CT(t)) Lt ) t — 0, 1, . . . . (10) 



As a consequence of the update based on the measure , the estimation 

-1 



error becomes = 



I-p-lCT 



It-l-l) 1 C'cr(t+l) 









c, 



a{t+l) 



-t+1 



In order to prove that the resulting estimator is stable, consider the Lyapunov 
functions Vt = ef P e^ and P , where P is a symmetric positive 

definite matrix. The goal is to demonstrate that the estimation error ef con- 
verges asymptotically to zero by proving that is decreasing in t, Ve+ G K" . 
To this end, it is sufficient to demonstrate that < Vt+i , t = 0, 1, . . . as, 
from m, it is obvious that Vj+i < , Ve+ G M” if the Lyapunov inequalities 

m are satisfied (see also the proof of Theorem Q]) . Thus, let us consider 






^t+l~P ^ (pcr{t+l) P C'cr(t-l-l) Lt-l-1 P 



-t+1 



-P 



-1 



P'^(t+l){po(t+l)P Cry(t+l)§Lt+l 



= (^(*+1) P-1 CT 



(i-l-l)) <^<7(4+1) 
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and we can conclude since <^^(4+1) (c'-r(t+i) P ^ C'J(t+i)) C^^t+i) ^t+i < 0 , 

Vet+i e K” . 



The design of the observer Q can be accomplished by solving the related 
Problem 0 in order to satisfy Q as the requirements to apply the projection 
update are only that P is a positive definite matrix and the matrices Ci G 
^ i = 1, 2, . . . , fc , are of full rank. 

The projection method has been successfully applied to the estimation 
of a class of continuous-time nonlinear systems with asynchronous measure- 
ments HH. Moreover, the performance improvements provided by the projection 
method will be highlighted by means of the simulation results presented in the 
next section. 



4 Simulation Results 

In order to show the effectiveness of the proposed estimation methods and the 
feasibility of the related LMI design procedure, let us consider the simple me- 
chanical system depicted in Fig. 01 



b, 




Fig. 4. Simple mechanical system. 



The continuous-time dynamics of the system is given by the following equa- 
tion: 

( i{t) = A x(t) + B u{t) , . 

\y(t) = C.p)x(t) 

where x{t) = [xi{t),X2{t),X3{t),X4{t)] G R"' is the state vector and u{t) G R 
is the scalar input. More specifically, Xi{t) is the position and of the mass mi , 
X2{t) is the speed of mi , x^it) is the position of the mass m2 , and 2:4 (t) is 
the speed of m2 . The matrices A and B are as follows: 





/ 0 1 0 1 \ 




( ^ \ 


A = 


— {ki + k2)lmi—bifmi—k2fmi 0 




1/mi 




0 0 0 1 




0 




y fc2/m2 0 — fc2/m2 — 62 /m 2 / 




1 0 / 
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where the parameters are k\ = 0.0642 , ^2 = 0.1925 , mi = 

5.1962 m2 = 1.7321 6i = 0.8660X5/3, and 62 = 0.1732 Xg/s. 
The system is switching in that the matrix C'g.(t) can assume values in the set 
{Ci,C2,C3,C4} with Cl = (1,0, 0,0,), C2 = (0,1, 0,0), C3 = (0,0, 1,0), and 
C4 = (0,0,0,l). The choice among the four candidate measurement equations is 
random, with the same probability of occurrence. The system is observable with 
any of the four measurement equations. The input u{t) was taken to be equal to 
a sinusoidal force, i.e., u{t) = ku sin (wt ) , where, for each simulation run, ku 
and w were randomly chosen in [0.0, 4.0] N and [0.1, 0.6] rad/s, respectively. 
Moreover, the initial states were randomly Gaussian distributed around 0 with 
standard deviations 5.0, 2.0, 5.0, and 2.0 for Xi,X 2 ,xs, and X 4 , respectively. 

A corresponding discrete-time model was obtained for the same system dis- 
cretizing equation (HU by means of a simple Euleur’s approximation with a time 
step equal to 0.1s. In the discrete-time setting, the standard routines of the 
Matlab LMI Control Toolbox provided the following solution to Problem 0 



/ 0.1331 0.1175 -0.0873 -0.0159 \ 
0.1175 1.4445 -0.0371 -0.1386 
-0.0873 -0.0371 0.0960 0.0351 

\ -0.0159 -0.1386 0.0351 0.8227/ 



/ 1.0028 \ 




/-1.5816\ 




/ 0.6846 \ 




/ 0.4669 \ 


-0.0623 




1.0246 




-0.0325 


u = 


0.1131 


0.9016 


L2 = 


-1.1822 


^3 = 


1.0065 


-0.6561 


\ -0.0280 ^ 




0.1874^ 




1^-0.0386 ^ 




0.9940/ 



The root mean square (RMS) error was considered as a performance index. This 
error for the scalar variable Xi(t) with respect to its estimate Xi(t) at the time 
t for N different trials is defined as 



RMSiit) = 



N 




xlit) - xlit) 

_____ 



(12) 



where x\ (t) is the value of the variable in the j-th run, xf (t) is the estimate of 
x/(t) , and i = 1,2,3, and 4. In Fig. El the simulation results obtained with the 
two proposed observers using the above-written gains and with initial estimated 
state equal to 0 are shown as far as it regards the RMS estimation error on 
500 trials (i.e., N = 500 ) with different choices of ku , w , initial state vectors, 
and switching sequences. As can be noticed in Fig. El the enhanced Luenberger 
observer exhibits a faster convergence rate. The trajectories of the true and 
estimated state variables for a single random-chosen simulation run are shown 
in Fig. El 



Position [m] Position [m] Position [m] Position [m] 
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RMS error on the estimate of x1 




RMS error on the estimate of x3 




Time [s] 



RMS error on the estimate of x2 




RMS error on the estimate of x4 




Fig. 5. RMS estimation errors of the switching observers. 



State variable x^ State variable x^ 





State variable x^ 




State variable x. 

4 




Time [s] 



Fig. 6. True values and estimates of the switching observers for a single realization. 
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5 Conclusions 

In this paper, estimation for a class of hybrid systems has been considered. First, 
we have addressed the problem of designing a Luenberger observer for a class of 
switching discrete-time linear systems. Conditions ensuring the stability of the 
error dynamics for such an estimator have been found and an LMI formulation 
has been presented to synthesize the gains in a straightforward, efficient way. 
Second, an enhanced Luenberger observer has been proposed to perform esti- 
mation using also the last available measures. The stability of the estimation 
error for this modified Luenberger observer has been proved under conditions 
that can be ensured by solving the same LMI problem of the first estimator. The 
simulation results obtained with such observers for a simple mechanical system 
show both that, as expected, the proposed estimators are stable and that the 
enhanced Luenberger observer results in higher performance. 

Future work will concern the application of the proposed approach to real 
cases (see M). where conventional estimation methods based on Kalman filter- 
ing may perform poorly. Moreover, further theoretical investigations will regard 
the extension of the switching observer to a more general framework, e.g., with 
noises acting on the system and measurement equations and nonlinearities af- 
fecting the dynamics. 
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Abstract. In a biological cell, cellular functions and the genetic regula- 
tory apparatus are implemented and controlled by a network of chemical 
reactions in which regulatory proteins can control genes that produce 
other regulators, which in turn control other genes. Further, the feed- 
back pathways appear to incorporate switches that result in changes in 
the dynamic behavior of the cell. This paper describes a hybrid systems 
approach to modeling the intra-cellular network using continuous differ- 
ential equations to model the feedback mechanisms and mode-switching 
to describe the changes in the underlying dynamics. We use two case 
studies to illustrate a modular approach to modeling such networks and 
describe the architectural and behavioral hierarchy in the underlying 
models. We describe these models using Charon P|, a language that 
allows formal description of hybrid systems. We provide preliminary sim- 
ulation results that demonstrate how our approach can help biologists 
in their analysis of noisy genetic circuits. Finally we describe our agenda 
for future work that includes the development of models and simulation 
for stochastic hybrid systemsQ 



1 Introduction 

In order to survive, organisms continuously monitor their surroundings and, if 
necessary, adjust traffic through simple or complex combinations of genetic and 
metabolic networks to respond to alterations in local conditions. Local condi- 
tions include both the physical environment, for example, temperature (the heat 
and cold shock response), nutrient and energy source concentrations (the strin- 
gent response), light (circadian rhythms), cell density (quorum sensing response) 
as well as the molecular environment of individual regulatory components. Ex- 
amples of the latter include intracellular concentrations of transcription factors 
and allosteric effectors. The availability of complete genomic information for a 
wide variety of organisms and the consequent attention on proteomics has dra- 
matically increased the number of systems and components of systems that are 
involved in these sensing and responding activities mni. Understanding how 

^ This research was supported in part by DARPA/ITO Mobies project (grant number 
F33615-00-C-1707) and NSF grant CDS-97-03220. 
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these biological systems are integrated and regulated and how the regulation 
may be influenced, possibly for therapeutic purposes, remains a significant chal- 
lenge. 

In this paper we model and simulate examples of genetic and metabolic net- 
works using a hybrid systems approach that combines concepts and tools from 
control theory and computer science. First we analyze a previously published 
plasmid-based genetic network that was designed and synthesized using three 
repressor transcription factors where one repressor negatively regulates the pro- 
duction of a subsequent repressor \J\. Then we model a biologically important 
genetic network that controls the quorum sensing response, an adaptive response 
of certain gram negative bacteria to local population density HM. The quorum 
sensing response controls the luminescent behavior in certain strains of Vibrio 
which has been linked to the normal development of the bacterial host HH| as 
well as to medically important phenomena such as biofllm formation by Pseu- 
domonas aerugenosa, an organism that can cause overwhelming pneumonia and 
septic shock mEDi. 

2 Modeling 

The genetic circuits and biomolecular networks considered here and elsewhere 
are remarkably similar to hybrid systems encountered in engineering, for exam- 
ple embedded systems. In particular, it is worth noting the following three key 
features: 

Concurrency and communication. At the intra-cellular level, proteins and 
mRNAs are agents communicating with each other and influencing each 
other’s behavior. At the inter-cellular level, cells can be viewed as networked 
agents interacting with each other via different communication mechanisms. 
Discrete and continuous behaviors. At the lowest level, the evolution of en- 
tities such as proteins can be described by differential equations. Discreteness 
arises in two ways. First, a certain activity may be triggered only when the 
concentration of enabling quantities is above the desired threshold. This leads 
to discrete switching between active and dormant states. Second, different 
models may be appropriate at different levels of concentration. 

Stochastic behavior. Evolution of entities is not deterministic, and is better 
captured by stochastic models that allow for uncertainty and noise. 

These characteristics are typical of high-level models of embedded software such 
as autonomous communicating mobile robots. For describing such systems, we 
have developed the language Charon |2| which incorporates ideas from con- 
currency theory (languages such as CSP m), object-oriented software design 
notations (such as Statecharts [3 and UML P|), and formal models for hybrid 
systems (such as hybrid automata P and hybrid I/O automata ^S|). The key 
features of Charon are: 

Architectural hierarchy. The building block for describing the system ar- 
chitecture is an agent that communicates with its environment via shared 
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variables. The language supports the operations of composition of agents to 
model concurrency, hiding of variables to restrict sharing of information, and 
instantiation of agents to support reuse. 

Behavior hierarchy. The building block for describing flow of control inside an 
atomic agent is a mode. A mode is basically a hierarchical state machine, that 
is, a mode can have submodes and transitions connecting them. Variables 
can be declared locally inside any mode with standard scoping rules for 
visibility. Modes can be connected to each other only via well-defined entry 
and exit points. We allow sharing of modes so that the same mode definition 
can be instantiated in multiple contexts. Finally, to support exceptions, the 
language allows group transitions from default exit points that are applicable 
to all enclosing modes. 

Discrete updates. Discrete updates are specified by guarded actions labeling 
transitions connecting the modes. Actions can have calls to externally defined 
Java functions which can be used to write complex data manipulations. It 
also allows us to mimic stochastic aspects through randomization. 
Continuous updates. Some of the variables in Charon can be declared ana- 
log, and they flow continuously during continuous updates that model pas- 
sage of time. The evolution of analog variables can be constrained in three 
ways: differential constraints (e.g. by equations such as i = f{x,u)), alge- 
braic constraints (e.g. by equations such as y = g{x,u)), and invariants (e.g. 
\x — y\ < s) which limit the allowed durations of flows. Such constraints can 
be declared at different levels of the mode hierarchy. 

Modular features of Charon allow succinct and structured description of 
complex systems. Similar features are supported by the languages Shift 0 and 
Stateflow (see www.mathworks.com). In Charon, modularity is not only ap- 
parent in syntax, but we are developing analysis tools (such as simulation) that 
exploit this modularity. Furthermore, Charon has formal foundations support- 
ing compositional refinement calculus which allows relating different models of 
the system in mathematically precise manner. A formal mathematical descrip- 
tion allows us to develop tools for computing equilibria, for reachability analysis 
and for analyzing properties like stability and reachability. 

In the next two sections, we will briefly describe case studies that we have 
used to investigate the hybrid systems approach to modeling biological systems, 
and the applications of Charon to these systems. We will also illustrate our 
approach by providing preliminary simulation results. 



3 A Repressilator Network 

As noted in most biomolecular systems of interest involve many interactions 
connected through positive and negative feedback loops and an understanding of 
their dynamics is hard to obtain. In this section we will describe the modeling of 
a specific biomolecular network. We will model a repressilator system described 
in First we provide some biological background information and describe the 
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protein network used in |3 , and then describe the models of the protein network, 
including examples of Charon modelsfl 

3.1 The Basic Phenomena 

In the synthetic oscillatory network described in [Z|, networks of interacting 
biomolecules carry out many essential functions in living cells. But the design 
principles of the functioning of such networks still remain poorly understood- 
even in relatively simple systems m. The authors proposed the design and 
construction of a synthetic protein network implementing a particular function. 
Their motivation is that such “rational network design” may lead to the engi- 
neering of new cellular behaviors and to improved understanding of naturally 
occuring networks. 

The repressilator system described in 0 contains three proteins, namely 
lad, tetR, and cl. The protein lad represses the protein tetR, tetR represses cl, 
whereas cl represses lad, thus completing a feedback system called a repressi- 
lator system. The dynamics of the network depend on the transcription rates, 
translation rates, and decay rates of proteins and messenger RNAs. Depending 
on the values of the different parameters in the model, the system might converge 
to a stable limit cycle or become unstable. 

3.2 Approaches to Modeling 

It is well known in mechanics and thermodynamics that there are two different 
approaches to modeling systems such as the repressilator system. At reasonably 
high molecular concentrations, one can adopt continuum models which lend 
themselves to deterministic models involving ordinary and partial differential 
equations. At lower concentrations, the discrete molecular interactions become 
important and deterministic models are difficult to obtain 0. 



The Deterministic, Continuous Approximation. We will consider the 
three repressor protein concentrations pi,i € P — {lad, tetR, cl} and their cor- 
responding mRNA concentrations G P as continuous dynamic variables. 
The system kinetics are determined by the following six coupled first-order dif- 
ferential equations. 



drui 

dt 

dpi 

dt 



-mi + — — - + oo 
1 +p” 



~(}{pi - m{) 



(i,j) G |(lacl,cl), (tetR,lacI), (cI,tetR)| 

^ For more information on Charon or sample Charon code, please check 
http://www.cis.npenn.edu/mobies/charon/ or contact ivancic@seas.upenn.edn. 
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The equations use various constants. The leakiness of the promoter a is the 
number of protein copies per cell produced from a given promoter type during 
continuous growth in the presence of saturating repressor amounts. During the 
absence of the repressor, we have a + ao number of protein copies per cell. The 
ratio of the protein decay rate to the mRNA decay rate is denoted by j3, while 
n stands for the so called Hill coefficient. 

The Stochastic, Discrete Approximation. The continuous analysis neglects 
the discrete nature of molecular components and the stochastic character of their 
interaction |Z]. Following [7], we adopt the stochastic approximation as described 
by Gillespie in |H|. The various proteins and mRNAs are modeled by discrete 
variables corresponding to the number of molecules measuring concentration, 
and are updated at discrete time intervals by stochastic rules. 

3.3 Charon Model 

In this section we will present the repressilator system models as described in | 7 ] 
using the Charon language. We will present many of the advantages that the 
Charon language has to offer for modeling such biomolecular models. 

Our model will define a generic protein model as an agent in Charon. We will 
instantiate this agent model to obtain the three proteins lad, tetR, and cl. The 
approximation models will be implemented inside the modes of the protein agent. 
To present another feature of our language, we will also describe a combination 
of the discrete and the continuous model into one modeling system. 

The Protein Agent in the Continuous Approximation. In this section 
we will describe a Charon model of a generic protein agent. We have a con- 
tinuous input variable which represents the repressor protein concentration pfi. 
This means, that the environment of this protein agent supplies the value of this 
variable, and it cannot be changed by the protein agent. The protein agent has a 
continuous private variable representing the messenger RNA concentration. Pri- 
vate variables cannot be seen outside the agent and can be updated internally for 
internal use only. The output of the protein agent is a continuous variable rep- 
resenting the protein concentration. Output variables are updated by the agent, 
and can be used as input variables to other agents. The generic protein agent 
has parameters ag, a, P,n,po, and mg. By instantiating these parameters with 
values, we can obtain instantiated protein agents representing a specific protein. 
The parameters po and ttiq will be used for initialization purposes and stand for 
the initial protein concentration and the initial messenger RNA concentration 
respectively. The following represents the corresponding Charon code. 

agent contProtein (real pO , mO , alphaO , alpha , beta , n){ 
write analog real p = pO ; //protein concentration 
read analog real pR ; //repressor protein concentration 
private ainalog real m = mO ; //messenger RNA concentration 
mode cent = continuous ( alphaO , alpha , beta , n ) ; } 
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read pR 



protein agent 

local m 

parameters; pO,mO,alphaO, alpha, beta, n 

P=pO , 
m=mO 



continuous mode 



d(m)=-m+alpha/( l+pR'^n)+alphaO 
d(p)=-beta*(p-m) 



Fig. 1. A generic protein agent for the continuous approximation model 



We still need to define the behavior of the agent. The behavior is described 
by the modes of the agent. The behavior of the generic protein agent is defined 
in cont, which is an instantiation of a generic continuous mode defined by the 
following code. A graphical version of the generic protein model can be found in 

Figure n 

mode continuous (real alphaO , alpha , beta , n){ 
write analog real p ; //protein concentration 
read analog real pR ; //repressor protein concentration 
private amalog real m ; //messenger RNA concentration 
diff mRNA { d(m) = -m + alpha / (l+pR"n) + alphaO } 
diff proteinConcentration { d(p) = -beta * (p-m) } } 




Fig. 2. Composed repressilator system using the instantiated generic protein agent 



Instantiation and Concurrency. We defined a generic protein agent in the 
previous section. We have to instantiate this generic agent model to get the 
three proteins used in the system. We also want the three proteins lad, tetR, 
and cl to run in parallel and to influence each other. Notice the use of renaming 
of variables to couple the three instantiated protein agents to influence each 
other. A graphical version of the composed system is illustrated in Figure|21 The 
following represents the corresponding Charon code using some values for the 
parameters. A simulation trace generated by the Charon tool-set is given in 
Figure 0 
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agent RepressilatorSystem (){ 



private ainalog real pi , 


P2 


p3 














agent lad = contProtein 


( . 


. ) 


[ 


P 


pR 


:= pi 


p3 ] 




agent tetR = contProtein 


( . 


. ) 


[ 


P 


pR 


:= p2 


pi ] 




agent cl = contProtein 


( . 


. ) 


[ 


P 


pR 


;= p3 


p2 ] 


} 




Fig. 3. Simulation trace for the repressilator system showing stable oscillations for the 
three protein concentration pi,P2,P3 over time. 



The Protein Agent in the Discrete Approximation. In this section we 
will present a possible model for a discrete approximation of a protein agent. As 
we did it for the continuous case, we will again define a generic protein agent, 
that can be instantiated to build a system of proteins. Our model works as 
follows. We have an integer variable n that keeps track of the number of protein 
molecules which is the output of the agent . The input to the agent is the number 
of repressor protein molecules ur. Depending on various parameters, we want 
to increase or decrease the number of protein molecules by one at a time. The 
basic idea is to use stochastic simulation as described in jS|. The parameters 
that influence the stochastic simulation are binding and unbinding of proteins 
on two-sided promoters, the protein and mRNA decay rates, and translation. 
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read pR 



local m 



) 1 






continuous mode 







discrete mode 




Fig. 4. A generic protein agent for the combined framework using continuous and 
discrete approximation model 



Combining the two Models into one Framework. The two different models 
for the repressilator system can be combined into one framework. The basic idea 
is to use the deterministic continuous model whenever the concentration of the 
protein is high enough, whereas we would switch to the discrete, stochastic model 
if the concentration would fall below a certain threshold value. Figure E| gives an 
intuitive graphical representation of the protein agent with both the continuous 
and discrete approximation. 



4 Quorum Sensing in Bacteria 

A good illustration of multicellular behavior in prokaryotes is the cell-density- 
dependent gene expression. In this process, a single cell is able to sense when 
a quorum of bacteria, a minimum population unit, is achieved. Under these 
conditions, certain behavior is efficiently performed by the quorum, such as bio- 
luminescence, which is the best known model for understanding the mechanism 
of cell-density-dependent gene expression. In this section, we will describe a hy- 
brid system model that captures the changes in dynamics of the biochemical 
reactions observed in the literature \rmsm- 



4.1 The Basic Phenomena 

Vibrio fischeri is a marine bacterium that can be found both as a free-living 
organism and as a symbiont of some marine fish and squid. As a free-living or- 
ganism, V. fisheri exists at low densities (less than 500 cells per ml of seawater) 
and appears to be non-luminescent. As a symbiont, the bacteria live at high 
densities and are, usually, luminescent. In a liquid culture, the bacteria’s level of 
luminescence is low until the culture reaches mid to late exponential phase. A 
dramatic increase in luminescence is observed at that time due to the transcrip- 
tional activation of the lux genes. Once the bacteria reach stationary phase, the 
level of luminescence decreases. 

The Zitxregulon im contains two operons, Ol and Or (see Figure EJ. The left 
operon Ol contains the luxR gene encoding the protein LuxR, a transcriptional 
activator of the system. The right operon Or contains seven genes luxICDABEG. 
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Fig. 5. A portion of DNA emphasizing luxR and luxICDABEG genes and the binding 
sites for LuxR complex and CRP 



Protein LuxI, the product of the luxi gene is required for endogenous production 
of autoinducer, a small molecule capable of diffusing in and out of the cell mem- 
brane. Genes luxA and luxB encode two subunits of luciferase. The trio luxC, 
luxD, and luxE code for the subunits of a protein complex which provides an 
aldehyde substrate for luciferase. The function of luxG is unknown. The autoin- 
ducer Ai binds to protein LuxR to form a complex Co. The two operons are 
separated by a regulatory region that contains a binding site for the cyclic AMP 
receptor protein CRP and a binding site for the complex Co. 

The transcription of luxR is regulated by both CRP and Co. We can distin- 
guish among the following three different cases: 

— Case Ol-1 In the absence of the autoinducer, CRP activates Ol expression 
by initiating two RNA transcripts. 

— Case Ol-“^ At low autoinducer concentrations, luxR transcription is stimu- 
lated by increasing CRP-dependent transcription and by Co-dependent tran- 
scription from another transcriptional start site. 

— Case Ol-3 At high autoinducer concentrations, luxR transcription is re- 
pressed through a second, weaker Co binding site located in luxD. 

Likewise, transcription of Or is regulated by both CRP and Co. We distinguish 
two different cases: 

— Case Or-1 In the absence of autoinducer, CRP represses Or transcription. 

— Case Or-2 In the presence of autoinducer, Co activates transcription of Or. 

These cases will be interpreted as modes as seen later in the paper. 

4.2 Mathematical Model 

In this section, we develop a mathematical model for the luminescence phe- 
nomenon in one bacterium of V. fischeri, describing the concentrations of the 
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relevant mRNA’s, proteins, and small molecules. As described in Section B~Tl the 
mechanism of transcription activation of both operons is highly dependent on 
the concentration of autoinducer, so the time evolution of the system cannot be 
described by one set of continuous differential equations^ Combining cases for 
Ol and Or given in the previous section, yields three modes, which we call OFF, 
POS and NEG. The transitions between modes are governed by the level of inter- 
nal autoinducer which we represent by [Ai]. Mode OFF corresponds to very low 
or zero concentration of autoinducer ([Ai] < [Ai]_) within the bacterium and no 
luminescence is observed. The system is in mode POS when the concentration 
of internal autoinducer is low ([Ai]_ < [Ai] < [Ai]_j_). This mode corresponds 
to positive growth and increasing concentration of autoinducer. Luminescence 
is observed, as are higher concentrations of proteins LuxA, LuxB, LuxC, LuxD, 
and LuxE. The transition to mode NEG (negative growth) occurs at high levels 
of autoinducer ([Ai] > [Ai]_j_). 

We use the following rate equation to describe the concentration for any 
molecular species (mRNA, protein, protein complex, or small molecule) m-- 



d\x^ 

= synthesis — decay ± transformation ± transport 



( 1 ) 



The synthesis term represents transcription for mRNA and translation for pro- 
teins. The decay term represents a first order degradation process. The transfor- 
mation term describes reactions such as cleavage or ligand-binding that do not 
destroy the protein, but do remove its ability to participate in specific reactions. 
Finally, molecular species may participate in transport processes, like passive 
diffusion or active transport through a membrane. 



The biomolecular system can be described in a nine dimensional state space. 
The nine variables, x\, X2, ■ ■ ■ , xg, describe the concentrations of different mole- 
cules as follows: 



Xi = mRNA transcribed from Or, 

X2 = mRNA transcribed from Or, 

X3 = protein LuxR, 

X4 = protein LuxI, 

X5 = protein LuxA, 

Xq = protein LuxB, 

X7 = autoinducer inside the bacterium Ai, 
xs = LuxR:Ai complex Co, 
xg = autoinducer outside the bacterium Aigx, 
where Ai is the dimensionless version of [Ai]. 



For simplicity, we have assumed that the concentrations of CRP and of the 
substrate necessary for endogenous production of Ai are constant. Further, we 
have neglected the decay rates for chemical compounds. Finally, we assume that 



^ In |1 H| . the differential equations for the low autoinducer concentration are described. 
The model presented here describes a wider range of operating conditions. 
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the concentrations of LuxC, LuxD, and LuxE are similar to those of Lux A and 
LuxB. 

The (continuous) differential equations for each mode are of the form 
X = p(x) where x = [xi, X2, ■■.,xg]'^ € IR®, P = [/(, P^, •.•,/g], and 
i G {OFF, POS, NEG{. The components of the vector fields are explicitly 
given by: 



fPOS 

Jl 



pNEG 

Jl 

pOFF 

/2 

pPOS _ pNEG 
J 2 — J2 



fl 

fl 

fl 

fl 

fl 



4 



3c - 



— 4a; 1 



^81 



-rjixi 

-V2X2 



m 



^82 

8 



^82 



+ X 



8 




773 (xi - 0:3) - r 3 T^AiX 3 .X 7 
r]4 {X2 — Xp — r^Xi 
V5 {X2 - X5) 

rje {X2 - xq) 



-f]7X7 + V4X4 - Tmem {x7 ~ Xg) - r37^RX3X7 



fl = -mxs + r37,AiX3X7 

fl = -V 7 X 9 + rmem{X 7 - Xg) + U 

where, in the last seven equations /* is independent of the mode. All the quan- 
tities in the above model are non-dimensional. r]i = Tg/Fli where Tg is the 
characteristic time constant of the system and Fli is the half-life (inverse of the 
decay rate) of molecule Xi- Vij is a cooperativity coefficient while describes 
the potency of the regulation of the transcription of mRNA j by protein i. r de- 
notes transformation and transfer rates. For example Tmem is the transfer rate of 
autoinducer through the membrane of the cell while r37^n and X37^Ai are transfor- 
mation rates obtained by non-dimensionalizing the binding rate of the reaction 
between Ai and LuxR in two different ways, c is dependent on the concentration 
of CRP and its affinity to the corresponding binding site, and, as stated ear- 
lier, is assumed to be constant. Finally, u emulates an external source of Ai and 
is used to simulate the sensitivity of the bacterium to changes of autoinducer 
concentration in the exterior. 

We regard u as an input to our system. Since proteins Lux A and LuxB are 
subunits of luciferase, which produces luminescence, it is reasonable to assume 
that the level of luminescence is proportional to the product of the concentrations 
of LuxA and LuxB, which we choose to be the output of the system. 



4.3 Charon Model 

The behavioral hierarchy in Charon (see Figure EJ is characterized by three 
different behaviors which are represented by three different modes, namely OFF, 
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POS, and NEG. Many of the differential equations governing the dynamics of 
the system are shared between the modes. We will introduce the notion of mode 
hierarchy to extract the shared constraints. Through the notion of submodes 
and scoping, we can simplify the description of the respective modes OFF, POS, 
and NEG. 




Fig. 6. Charon structure of the system 



Figure 13 illustrates the response {i.e., luminescence) of the bacterium to a 
perturbation in the concentration of external autoinducer that takes the form 
of a rectangular pulse. The magnitude of the step has been chosen to make 
the system go through all three modes. The results confirm the experimental 
observations m luminescence increases during mode POS and decreases in 
mode NEG; there is no luminescence in mode OFF. The switch history and the 
time evolution of the concentrations of the significant molecules in the system 
are also shown. In mode OFF, all molecules decay to zero, except for mRNA 
Ol and the corresponding protein R, as expected. For a short time, in mode 
POS, all the concentrations increase until the internal autoinducer reaches a high 
concentration, when the system is switched to mode NEG. In this last mode, 
everything decays to zero, except for internal autoinducer which can reach a 
stable non-zero value dependent on the size of the step of external autoinducer. 



5 Conclusions 

In this paper we have shown that biological cellular networks can be natu- 
rally modeled as hybrid systems. In particular, the protein repressilator system 
switches between a continuous deterministic model at high concentrations, and 
a timed, discrete, stochastic model at low concentrations. Similarly, the lumi- 
nescence control of Vibrio fischeri is naturally modeled as a multi-modal hybrid 
system, resulting in simulations that are in accordance with experimental obser- 
vations. The hybrid nature of such protein networks can be very easily expressed 
and simulated in Charon, which may offer us better and a more global under- 
standing of biological networks. 
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Fig. 7. Increase in external autoinducer produces luminescence: (a)input - external 
source of autoinducer; (b) switch history; (c) output (luminescence)- product of con- 
centrations of proteins A and B; (d) and (e) time - evolution of concentrations; 



The enormous complexity of large scale biological networks will present us 
with great challenges that we must face. Exploiting the structure of biological 
systems will be critical for scaling the applicability of the modeling, analysis, and 
simulation tools. It is therefore extremely encouraging that the two case studies 
presented in this paper exhibit the architectural paradigms of modern software 
engineering. 

We envision the link between hybrid systems technology, and biology to 
strengthen. The scalable nature of computational tools like Charon will en- 
able the unified and improved modeling of biological cellular networks, leading 
to better understanding, as well as providing us with the opportunity to deter- 
mine how local biological changes can affect global behavior. Conversely, a good 
understanding of the robustness of noisy biological networks will lead to new 
approaches to designing networked embedded systems. 

The case studies also highlight the need for developing a theory of stochastic 
hybrid systems, for instance, for modeling rate equations of biochemical reac- 
tions. We believe that mathematical and computational tools for the analysis 
of such systems present a research challenge for the hybrid systems commu- 
nity, while presenting a significant potential for greatly impacting post genomics 
research. 



References 

1. R. Alur, C. Courcoubetis, N. Halbwachs, T.A. Henzinger, P. Ho, X. Nicollin, A. 
Olivero, J. Sifakis, and S. Yovine. The algorithmic analysis of hybrid systems. 
Theoretical Computer Science, 138:3-34, 1995. 








32 



R. Alur et al. 



2. R. Alur, R. Grosu, Y. Hur, V. Kumar, and I. Lee. Modular specifications of hybrid 
systems in Charon In Hybrid Systems: Computation and Control, Third Interna- 
tional Workshop, volume LNCS 1790, pages 6-19, 2000. 

3. G. Booch, I. Jacobson, and J. Rumbaugh. Unified Modeling Language User Guide. 
Addison Wesley, 1997. 

4. R. Brent. Genomic biology. Cell, 100(1):169-183, January 2000. 

5. H. de Jong, M. Page, C. Hernandez, H. Geiselmann, and S. Maza. Modeling and 
Simulation of Genetic Regulatory Networks. ERCIM News, 43, October 2000. 

6. A. Deshpande, A. Gollu, and L. Semenzato. SHIFT programming language and 
run-time systems for dynamic networks of hybrid automata. Technical report. Uni- 
versity of California at Berkeley, 1997. 

7. M. Elowitz and S. Leibler. Asynthetic oscillatory network of transciptional regula- 
tors. Nature, 403:335-338, January 2000. 

8. D.T. Gillespie. Exact stochastic simulation of coupled chemical reactions. J. Phys. 
Chem., 81:2340-2361, 1977. 

9. D. Harel. Statecharts: A visual formalism for complex systems. Science of Com- 
puter Programming, 8:231-274, 1987. 

10. L.H. Hartwell, J.J. Hopfield, S. Leibler, and A.W. Murray. From molecular to 
modular cell biology. Nature, 402((6761 Suppl)):C47-52, December 1999. 

11. D.J. Hassett, J.F. Ma, J.G. Elkins, T.R. McDermott, U.A. Ochsner, S.E. West, 
C.T. Huand, J. Fredericks, S. Burnett, P.S. Stewart, G. McFeters, L. Passador, 
and B.H. Iglewski. Quorum sensing in Pseudomonas aeruginosa controls expression 
of catalase and superoxide dismutase genes and mediates biofilm susceptibility to 
hydrogen peroxide. Mol Microbiol, 34(5): 1082-1093, December 1999. 

12. C.A.R. Hoare. Communicating Sequential Processes. Prentice-Hall, 1985. 

13. S. James, P. Nilson, G. James, S. Kjellenberg, and T. Fagerstrom. Luminescence 
control in the marine bacterium Vibrio fischeri: An analysis of the dynamics of lux 
regulation. J Mol Biol, 296(4):1127-1137, March 2000. 

14. D.E. Jr Koshland. The era of pathway quantification. Science, 280:852-853, 1998. 

15. N. Lynch, R. Segala, F. Vaandrager, and H. Weinberg. Hybrid I/O automata. In 
Hybrid Systems HI: Verification and Control, LNGS 1066, pages 496-510, 1996. 

16. H. H. McAdams and A. Arkin. Simulation of prokaryotic genetic circuits. Annu. 
Rev. Biophys. Biomol. Struct., 27:199-224, 1998. 

17. D.M. Sitnikov, J.B. Schineller and T.O. Baldwin. Transcriptional regulation of bi- 
oluminesence genes from Vibrio fischeri. Mol Microbiol, 17(5):801-812, September 
1995. 

18. K.L. Visick, J. Foster, J. Doino, M. McFall-Ngai and E.G. Ruby. Vibrio hscheri 
lux genes play an important role in colonization and development of the host light 
organ. Bacteriol, 182(16):4578-4586, August 2000. 

19. G. von Dassow, E. Meir, E. M. Munro, and G. M. Odell. The segment polarity 
network is a robust development module. Nature, 406:188-192, July 2000. 

20. H. Yang, M. Matewish, I. Loubens, D.G. Storey, J.S. Lam, and S. Jin. migA, a 
quorum-responsive gene of Pseudomonas aeruginosa, is highly expressed in the cys- 
tic fibrosis lung environment and modifies low-molecular-mass lipopolysaccharide. 
Microbiology, 146((Pt 10)):2509-2519, October 2000. 




Compositional Refinement for Hierarchical 
Hybrid Systems’^ 



Rajeev Alur^, Radu Grosu^, Insup Lee^, and Oleg Sokolsky^ 

^ Department of Computer and Information Science, University of Pennsylvania 
^ Department of Computer Science, State University of New York at Stony Brook 



Abstract. In this paper, we develop a theory of modular design and 
refinement of hierarchical hybrid systems. In particular, we present com- 
positional trace-based semantics for the language Charon that allows 
modular specification of interacting hybrid systems. For hierarchical de- 
scription of the system architecture, Charon supports building complex 
agents via the operations of instantiation, hiding, and parallel composi- 
tion. For hierarchical description of the behavior of atomic components, 
Charon supports building complex modes via the operations of instan- 
tiation, scoping, and encapsulation. We develop an observational trace 
semantics for agents as well as for modes, and define a notion of refine- 
ment for both, based on trace inclusion. We show this semantics to be 
compositional with respect to the constructs in the language. 



1 Introduction 

Modern software design paradigms promote hierarchy as one of the key con- 
structs for structuring complex specifications. We are concerned with two dis- 
tinct notions of hierarchy. In architectural hierarchy, a system with a collection 
of communicating agents is constructed by parallel composition of atomic agents, 
and in behavioral hierarchy, the behavior of an individual agent is described by 
hierarchical sequential composition. The former hierarchy is present in almost all 
concurrency formalisms, and the latter, while present in all block-structured pro- 
gramming languages, was introduced for state-machine-based modeling in Stat- 
ECHARTS PI, and forms an integral part of modern notations such as UML p). 

A hybrid system typically consists of a collection of digital programs that 
interact with each other and with an analog environment. Specifications of hybrid 
systems integrate state-machine models of discrete behavior with differential 
equations for continuous behavior. This paper is about developing a formal and 
compositional semantics of hierarchical hybrid specifications. Formal semantics 
leads to definitions of semantic equivalence (or refinement) of specifications based 
on their observable behaviors, and compositionality means that semantics of a 
component can be constructed from the semantics of its subcomponents. Such 
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formal compositional semantics is a cornerstone of concurrency frameworks such 
as CSP ^ and CCS fT^ , and is a prerequisite for developing modular reasoning 
principles such as compositional model checking and systematic design principles 
such as stepwise refinement. 

The main contribution of the paper is a formal compositional semantics for 
the language Charon P] with an accompanying compositional refinement cal- 
culus. The building block for describing the system architecture is an agent that 
communicates with its environment via shared variables. The language supports 
the operations of composition of agents to model concurrency, hiding of variables 
to restrict sharing of information, and instantiation of agents to support reuse. 
The building block for describing flow of control inside an atomic agent is a 
mode. A mode is basically a hierarchical state machine, that is, a mode can have 
submodes and transitions connecting them. Variables can be declared locally in- 
side any mode with standard scoping rules for visibility. Modes can be connected 
to each other only via well-defined entry and exit points. We allow sharing of 
modes so that the same mode definition can be instantiated in multiple con- 
texts. To support exceptions, the language allows group transitions from default 
exit points that are applicable to all enclosing modes, and to support history 
retention, the language allows default entry transitions that restore the local 
state within a mode from the most recent exit. Discrete updates are specified by 
guarded actions labeling transitions connecting the modes. Some of the variables 
in Charon can be declared analog, and they flow continuously during continu- 
ous updates that model passage of time. The evolution of analog variables can 
be constrained in three ways: differential constraints (e.g. by equations such as 
X = f{x,u)), algebraic constraints (e.g. by equations such as y = g{x,u)), and 
invariants (e.g. \x — y\ < e) which limit the allowed durations of flows. Such 
constraints can be declared at different levels of the mode hierarchy. 

To define the modular semantics for modes, with each mode we associate two 
relations, one capturing its discrete behavior and one capturing its continuous 
behavior. Defining the discrete relation is tricky in presence of features such 
as group transitions, exceptions, and history retention. Our solution relies on a 
closure construction, inspired by a similar construction for hierarchical discrete 
systems |2| , which allows us to treat the transfer of control between a mode and 
its environment as a game. 

While discrete steps of a mode and its environment are interleaved, continu- 
ous steps need to be synchronized as time is a global parameter. In fact, during 
a flow, all active hierarchically nested modes must participate. To allow flexible 
and hierarchical specifications, in Charon, flow constraints can be specified at 
all levels of the hierarchy. To formalize this feature in a consistent and modular 
manner, we require that a mode can participate in a flow only when the control 
is at its default exit point. Then, all applicable constraints are properly used to 
define permitted flows. 

The discrete and continuous relations of a mode allow us to define executions 
of a mode, and corresponding traces are obtained by projecting out the private 
variables. We show that the set of traces of a mode can be constructed from 
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the traces of its submodes. This compositionality result leads to a compositional 
notion of refinement for modes. A mode M refines a mode N if they have 
the same interface in terms of entry/exit points and shared variables, and the 
traces of M is a subset of traces of N. This notion admits modular reasoning 
in the following manner. Suppose we obtain an implementation design I from 
a specification design S simply by locally replacing some submode in S' by a 
submode M. Then, to show I refines S, it suffices to show that M refines N. 
We illustrate this benefit by a simple example. 

Once we have the compositionality results for modes, analogous results for 
agents are relatively straightforward. We define an observational trace semantics 
for agents, a resulting notion of refinement, and show it to be compositional with 
respect to the operations of parallel composition, hiding, and instantiation. 

Related work. Early formal models for hybrid systems include phase tran- 
sition systems m and hybrid automata P- Models such as hybrid I/O au- 
tomata H2I and hybrid modules ^ allow compositional treatment of concurrent 
hybrid behaviors. The notion of hierarchical state machines was introduced in 
Statecharts 0 , and is present in many software design paradigms such as 
Uml I^. Our treatment of hierarchy is closest to hierarchical reactive mod- 
ules PI which shows how to define a modular semantics for hierarchical (dis- 
crete) modes. Tools such as Shift [Zj, Ptolemy jS], and Stateflow (see 
www.mathworks.com) allow hierarchical specifications of hybrid behavior, but 
formal semantics has not been a concern. HyCharts |H| presents a hierarchical 
model with modular operational semantics, but does not consider refinement. 
Masaccio m is a formal model for hierarchical hybrid systems. While same in 
spirit, it differs from our model in many technically significant aspects: it allows 
nesting of sequential and parallel composition, and allows a more general form of 
synchronous communication, but disallows high-level features of Charon modes 
such as exceptions, history retention, and specification of constraints at various 
levels. 

2 Motivational Example 

In this section, we present a simple example that outlines features, useful in 
a specification language for hybrid systems. We also point out the difficulties 
of defining semantics for such a language. Then we give the intuition for our 
approach to the semantics definition, which allows us to overcome the difficulties. 

Our example is a system that controls the level of liquid in a leaky tank. 
The level is controlled by infusing a flow of liquid into the tank. The level in the 
tank can be measured directly, but the rate of the leak has to be estimated. The 
controller has two goals: first, it must make sure that the level is within some 
critical bounds. If it is not, emergency measures are taken to make the level safe. 
When the level is safe, the controller should change the infusion rate according 
to instructions of the user. To do that, the controller periodically recomputes 
the desired rate of change for infusion and maintains the computed rate until 
the next update. 
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We now present a hierarchical description of the system in Charon. The 
hierarchy in Charon is twofold. The architectural hierarchy describes how the 
system agents interact with each other, hiding the details of interaction between 
sub-agents. The behavioral hierarchy describes behavior of each agent, hiding the 
low-level behavioral details. In our example, we have only one level of architecture 
description with agents Tank and Controller. There are two variables shared 
by the agents: level for the level of the liquid, and infusion for the infusion 
rate. 

Both agents are primitive, that is, without concurrent sub-agents. Behavior 
of a primitive agent is given by a mode, a hybrid state machine equipped with 
analog and discrete variables. While a mode stays in a state, its analog variables 
are updated continuously according to a set of constraints. Taking transitions 
from one state to another, the mode updates its discrete variables. States of the 
mode are submodes that can have their own behavior. A mode has a number 
of control points, through which control enters and exits the mode. That is, to 
perform a computation in one of its submodes, a mode takes a transition to an 
entry point of that submode. When the computation is complete, a transition 
from an exit point of the submode is taken. Before the computation of a mode 
is completed, it may be interrupted by a group transition, originating from a 
default exit point dx. After an interrupt, control is restored to the mode via a 
default entry point de. In our example, the behavior of Tank is represented by 
a single differential equation d{level) = infusion — leak, where leak is a local 
variable of Tank. Figure ^ shows the behavior of the agent Controller. The 
top-level mode of Controller has two submodes, Normal and Emergency. We 
do not show the details of the mode Emergency. It is activated when the level 
enters the critical region. 




Compute 

local discrete real esr 

global analog real level, infusion 

global discrete real rate 




ComputeHigh 

d(infusion) = est-1 



ComputeLow 
d(infiision) = est+I 




Fig. 1. Behavior of the controller 



The mode Normal has two submodes. Submode Maintain is used to maintain 
the current rate of change for infusion, represented by a local variable rate. Every 
10 seconds, measured by a local clock t. Maintain makes a call to Submode 
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Compute that computes a new value of rate. The details of the computation are 
irrelevant, but we assume that the computation is done differently depending on 
the level. We therefore introduce two submodes in Compute and show only the 
constraints for infusion in each submode. The exit transition of Compute assigns 
the computed value to the variable rate. 

Note that the mode Normal controls the value of the clock t, and its rate 
of change is the same in all its submodes. By contrast, infusion is updated 
differently in the two submodes. In this case, every submode must provide a 
constraint for infusion. Note also that rate is a discrete variable. It is updated 
only by transitions of Compute. 

We use invariants to force one of the outgoing transitions. Control can reside 
in a mode only as long as its invariant is satisfied. As soon as an invariant is 
violated, control has to leave the mode by taking one of the enabled outgoing 
transitions. In Figured, invariants of the modes are shown in braces. For exam- 
ple, ten time units after entering the mode Maintain the transition to Compute 
has to be taken. 

We distinguish between regular transitions and interrupts. For example, con- 
trol is transferred from Compute to Maintain only when the computation is com- 
plete. When it is time to perform another computation, it will start from the 
beginning. On the other hand, the transition from Normal to Emergency works 
as an interrupt. Regardless of which submode of Normal is operating when an 
interrupt occurs, control is transferred to Emergency. Upon return from the in- 
terrupt, the control state of Normal is restored. There is no priority between 
regular transitions and interrupt^. A mode can ignore an enabled interrupt and 
execute its internal transitions or let time elapse. We use invariants as described 
above to enforce interrupts (see the invariant of mode Normal). Invariants give 
the user finer control over interrupts. For example, a situation when an interrupt 
is optional for some time and then becomes urgent can be easily expressed. 

In addition to discrete steps described above, a mode can make continuous 
steps, when time progresses and the analog variables of the mode are updated 
according to a set of constraints. Because of the hierarchical structure of the 
mode, the set of applicable constraints consists of the constraints defined in the 
mode itself and those from the currently active submode. This implies that a 
mode can engage in a continuous step only when its control properly resides 
within one of its submodes. For example, we cannot allow time to pass at the 
control point e of Compute, between executing the transition from Maintain to 
Compute and a transition to enter ComputeHigh or ComputeLow. 



3 Modes 

Notation. We will represent modes and agents as tuples of components. If T is 
a tuple (ti, . . . ,t„), we identify the component ti of T as T.ti. We extend this 



^ Other treatments of interrupts can be handled equally well within the proposed 
framework. For example, [2| discuss weak interrupts in a similar setting. 
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notation to sets of tuples. If ST is a set of tuples with the same structure, we 
write ST.ti to mean T-tj. 

Given a set V of typed variables, a valuation for R is a function mapping 
variables to their values. We will assume that all valuations are type correct. 
The set of valuations over V is denoted Qv- We will use variables s,t, possibly 
primed or subscripted, to range over valuations. Given a valuation s over V, and 
a set W C V, s[W] denotes the restriction of s to the variables of W. 

A flow for a set V of variables is a differentiable function / from a closed 
interval of non-negative reals [0,5] to Qv- We refer to 6 as the duration of the 
flow. We assume that only constant functions are differentiable for non real- 
valued types. We denote a set of flows for V as Ty- 



3.1 Syntax 

Definition 1. (Mode) A mode M is a tuple {E,X,V, SM, Cons,T), where E 
is a set of entry control points, X is a set of exit control points, V is a set of 
variables, SM is a set of submodes, Cons is a set of constraints, and T is a set 
of transitions. 

Variables. A mode has a finite set of typed variables V, partitioned into subsets 
Va and Vd, the sets of analog and discrete variables, respectively. We also parition 
V into Vg and Vi, the sets of global and local variable^ We assume that there 
are no conflicts between the names of local variables of different modes. 
Submodes. SM is a finite set of submodes. We require that each global variable 
of a submode is a variable (either global or local) of its parent mode. That is, if 
N G SM, then N.Vg C V. This induces a natural scoping rule for variables in a 
hierarchy of modes: a variable introduced as local in a mode is accessible in all 
its submodes but not in any other mode. 

Control points. E is the set of entry points; X is the set of exit points. There 
are two distinguished control points representing default entry and exit: de G E 
and dx G X. We use C for the set of all control points of the mode: C = 
A U A U SM.E U SM.X. 

Constraints. The finite set Cons of constraints defines the flows permitted by 
Af0. Cons contains an invariant I, which defines when the mode can be active 
(see the definition of an active mode below). Further, for a variable x G Va, 
Cons can contain an algebraic constraint Aj,, which defines the set of admissible 
values for x, or a differential constraint Da,, which defines admissible values 
for the derivative of x with respect to time. Every invariant and an algebraic 
constraint is a predicate c C Qy and a differential constraint Dj, is a predicate 
on Qvud(v). A flow / is permitted by the mode if for every t in the domain of /, 
every variable in f(t) satisfies all constraints in Cons. Examples of constraints 
are d{x) < f{x,y) and g{x,y) < 0. 

^ Charon refines the set of global variables further according to allowed read/write 
access, but we won’t make such a distinction in this paper for clarity of presentation. 
^ The semantics does not depend on how sets of flows are specified. Here, we chose 
one of the possible ways. 
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Transitions. T is a finite set of transitions of the form (e,a,a;), where e G 
E U SM.X, cc G X U SM.E, and a, the action of the transition, is a relation 
from Qy^ to Qy if e G if and from Qy to Qy otherwise. A transition connects 
control points of the mode or its submodes. When a transition is executed, it 
updates some variables of the mode. Every mode is assumed to have an identity 
transition from de to dx, but we disallow transitions from any non-default control 
point to dx. A transition that originates at a default exit point of a submode is 
called a group transition of that submode. A group transition can be executed 
to interrupt the execution of the submode. We require that if a submode has 
been exited by a group transition, it must be entered again through its default 
entry point to resume the interrupted execution. 

Furthermore, we require that the mode cannot be blocked at any of its non- 
default control points. Precisely, for every e of M that is not de in M or dx in one 
of the submodes of M, the union Oe of all actions of the transitions originating 
at e is complete, that is, for every s there is t such that (s,t) G ae- 
Special modes. We distinguish two kinds of modes that play a special role in 
the semantic definitions. A mode M is a leaf mode if M.SM = 0. Leaf modes 
perform continuous steps according to their constraints. A top-level mode has 
a single non-default entry point init and no non-default exit points. Top-level 
modes are used to describe behavior of agents, as shown in Section 0 



3.2 Semantics 

Intuition. A mode can engage in a discrete or continuous behavior. During 
an execution, the mode and its environment either take turns making discrete 
steps or take a continuous step together. Discrete and continuous steps of the 
mode alternate. During a continuous step, the mode follows a flow from the set 
of flows possible for the current state for the length of its duration, updating 
its variables according to the flow. Note that the set of flows permitted by the 
mode’s constraints may be further restricted by the mode’s environment. A 
discrete step of the mode is a finite sequence of discrete steps of the submodes 
and enabled transitions of the mode itself. A discrete step begins in the current 
state of the mode and ends when it reaches an exit point or when the mode 
decides to yield control to the environment and let it make the choice of the 
next step. Note that in the latter case, the decision to break a discrete step is 
made by the mode itself. Technically, when the mode ends its discrete step in 
one of its submodes, it returns control to the environment via its default exit 
point. The closure construction, described below, ensures that the mode can 
yield control at appropriate moments, and that the discrete control state of the 
mode is restored when the environment schedules the next discrete step. 

State of a mode. We define the state of a mode in terms of all variables of the 
mode and its submodes. We use Vt = V L) SM.V^ for the set of all variables. 

The state of a mode M is a pair (c, s), where c is the location of discrete 
control in the mode and s G Qm.v,- Whenever the mode has control, it resides 
in one of its control points. In this case, c G M.C. We use special symbol e to 
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denote the case when the mode does not have control. Given a state (c, s) of M, 
we refer to c as the control state of M and to s as the data state of M. 
Preemption. An execution of a mode can be preempted by a group transition. 
A group transition of a mode originates at the default exit of the mode. During 
any discrete step of the mode, control can be transferred to the default exit and 
an enabled group transition can be selected. There is no priority between the 
transitions of a mode and its group transitions. When an execution of a mode is 
preempted, the control state of the mode is recorded in a special history variable, 
a new local variable that we introduce into every mode. Then, when the mode is 
entered through the default entry point next time, the control state of the mode 
is restored according to the history variable. 

The history variable and active submodes. In order to record the location 
of discrete control during executions, we introduce a new local variable h into 
each mode that has submodes. The history variable ft, of a mode M can assume 
values from the set SMUe. A submode A^ of M is called active when the history 
variable of M has the value N. Every top-level mode is always active. 

Closure of a mode. Closure construction is a technical device to allow the mode 
to interrupt its execution, either to allow the environment to schedule another 
step or to provide for preemption of the mode execution by group transitions. 
Transitions of the mode are modified to update ft after a transition is executed. 
In addition, default entry and exit transitions are added to the set of transi- 
tions of the mode. These default transitions do not affect the history variable 
and allow us to interrupt an execution and then resume it later from the same 
point. 

The closure modifies the transitions of M in such a way that, after each 
transition, ft records the active submode. If a transition leads to a control point 
of a submode N, the resulting state has h = N. Otherwise, if the transition 
leads to a control point of M itself, the value of ft after the transition will be 
e. For each submode N oi M, the closure adds a default exit transition from 
N.dx to M.dx. This transition does not change any variables of the mode and 
is always enabled. Default entry transitions are used to restore the local control 
state of M. A default entry transition leads from a default entry of the mode to 
the default entry of every submode N and is enabled if h = N. Furthermore, we 
make sure that the default entry transitions do not interfere with regular entry 
transitions originating from de. The closure changes each such transition so that 
it is enabled only if ft = e. 

Formally, the closure c(M) of a mode M = {E, X, V, SM, Cons, T) is defined 
to be the mode {E, X, FUft, c{SM), Cons, c{T)), where ft ^ F is a new local vari- 
able, c{SM) = {c(m) I m G SM} is the set of closed submodes of M, and c(T) is 
the closed set of transitions obtained by extending T with transitions {x, ax,dx) 
for every x € SM.dx and {de,ax,e) for every e G SM.de, and extending every 
transition in T such that 

~ (s, s) G ax iff X G N.E for some N G SM and s[ft] = N-, 

— for every transition (e, a, x) G T, the respective closed transition is (e, a' ,x), 
where (s,t) G a' iff {s\V],t\V]) G a and 
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— if X G N.E for some N G SM, then t[h] = N, otherwise t[h] = e, 

— if e G N.X for some N G SM, then s[/i] = N, otherwise s[h] = e. 

The closure construction for the example introduced in Section |^is illustrated 
in Figure 13 To avoid cluttering the figure, we omit the default transitions of the 
submode ComputeLow, and do not show the variables of the modes. 




Compute 
h=ComputeHigh 

ComputeHigh 
d( infusion) = est-1 
ComputeHigh 




rate - 

level > 5 A h = z e 

level < 5 A h = E 

rate - 

ComputeLow 



ComputeLow 
d( infusion) = est+1 



Fig. 2. Closed modes 



Before formally defining executions of a mode, we illustrate continuous and 
discrete steps using the example in Figure 0 Assume that the the controller 
is in the Maintain mode and none of the invariants is violated. Maintain can 
voluntarily relinquish control to the environment to let it take a step or advance 
time by taking the default exit transition to dx of Normal. There, the group 
transition is not enabled, and the default exit transition of the parent mode 
is taken. When the control arrives thus at the top level, the environment can 
schedule a continuous step. The analog variables of all agents are updated ac- 
cording to the constraints of the active modes. The active modes are Maintain, 
Normal, and Controller. Thus, the applicable constraints are d{t) = 1 and 
d{in fusion) = rate. The global variable level is updated according to the con- 
straint in Tank. After the continuous step, control returns to Maintain via the 
chain of default entry transitions. Assume now that the invariant of Normal is vi- 
olated while control is inside a submode of Compute. Then, control is transferred 
to dx of Compute and then on to dx of Normal. There, the choice between the 
group transition to Emergency or the default exit transition is non-deterministic. 
But since the invariant is violated, a continuous step cannot be taken. 
Operational semantics. An operational view of a closed mode M with the set 
of variables V consists of a continuous relation and, for each pair c\ G E, 
C 2 G X, a discrete relation Rf)^c 2 - 

The relation R^ C Qy x Ey gives, for every data state of the mode, the set of 
flows from this state. By definition, if the control state of the mode is not at dx, 
the set of flows for the state is empty. We require that, whenever (s, /) G R^ , 
/(O) = s. In addition, for each s, the set of flows Eg = {f \ (s, f) G R^} is 
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prefix-closed. That is, if the domain of f € Ts is [0,(5], then for every e < 5, a 
flow /' : [0,e] that coincides with / on [0,e] also belongs to iFs- is obtained 
from the constraints of a mode and relations SM.R'^ of its submodes. Given a 
data state s of a mode M, (s, /) G R^ iff f is permitted by M and, if N is the 
active submode at s, {s[N.V], f[N.V]) G N.R^ . 

For each ci G if U SM.X, C 2 G X U SM.E, relation i?^_c 2 — Qv x Qv 
describes the discrete behavior in which control is transferred from ci to C 2 - 
The relation R^^ comprises macro-steps of a mode starting at e and ending 
at X. A macro step consists of a sequence of micro-steps. Each micro-step is 
either a transition of the mode or a macro-step of one of its submodes. Given 
the relations R^/ e' G SM.E, x' G SM.X of macro-steps of the submodes of 
M, a micro-execution of a mode M = {E, X, V, SM, C, T) is a sequence of the 
form (eo. So)) (ei, si), . . . ,{en,s„) such that, for all i, Ci G C and Si G 14 and 
for even i, ((cj, s^), (e^+i, s*+i)) G T, while for odd i, (si,Sj+i) G SM.R^. ,,.^^. 
Given such a micro execution of M with cq = e G E and e„ = a: G A, we have 
(so, Sji) G Re^x- 

Definition 2. (Operational semantics) The operational semantics of the mode 
M consists of its control points E 0 X, its variables V and relations R'^ and 

tdD 
^e,x ' 

The operational semantics of a mode defines a transition system TZ over 
the states of the mode. We write (ei, Si)-^(c 2 , S 2 ) if (si,S 2 ) G R^, b 2 ’ 

{dx, Si)^{dx, S 2 ) if (si,/) G R'^ , f is defined on the interval [0, t] and f{t) = S 2 . 
We extend TZ to include environment steps. An environment step begins at an 
exit point of the mode and ends at an entry point. It represents changes to the 
global variables of the mode by other components while the mode is inactive. 
Private variables of the mode are unaffected by environment steps. Thus there 
is an environment step (x, s)A(e,t) whenever x G X, e G E, and s[V^j = t\Vp\. 
We let A range over Ey U {o, e}. An execution of a mode is now a path through 
the graph of TZ: 



(eo) So)— )(ei, si)— I- . . . ^(e„, Sn). 



3.3 Trace Semantics 

To be able to define a refinement relation between modes, we consider a trace 
semantics for modes. A trace of the mode is a projection of its execution onto 
the global variables of the mode. That is, a trace is obtained from each execution 
by replacing every Si with sjVg], and every / in transition labels with f[Vg]. We 
denote the set of traces of a mode M by Lm . 

Definition 3. (Trace semantics for modes) The trace semantics for M is given 
by its control points E and X, its global variables V, and its set of its traces Lm. 

In defining compositional and hierarchical semantics, one has to decide, what 
details of the behavior of lower-level components are observable at higher levels. 
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In our approach, the effect of a descrete step that updates only local variables 
of a mode is not observable by its environment, but stoppage of time introduced 
by such step is observable. For example, consider two systems, one of which is 
always idle, while the other updates a local variable every second. These two 
systems are different, since the second one does not have flows more than one 
second long. Defining a modular semantics in a way that such distinction is not 
made seems much more difficult. 



4 Agents 

4.1 Syntax 

Definition 4. (Agent) An agent (TM,V,I) consists of a set of variables V, a 
set of initial states, and a set of top-level modes TM. 

The top-level modes collectively define behavior of the agent. The set V 
is partitioned into local variables Vi and global variables Vg. We require that 
TM.V C y, Rg C TM.Vg] that is, all global variables originate in some mode. 
The set of initial states I C Qy specifies possible initializations of the variables 
of the agent. A primitive agent has a single top-level mode. Composite agents 
have many top-level modes and are constructed by parallel composition of other 
agents as described below. 

4.2 Semantics 

An execution of an agent follows a trajectory, which starts in one of the initial 
states and is a sequence of flows interleaved with discrete updates to the variables 
of the agent. An execution of A is constructed from the relations and of 
its top-level modes. For a fixed initial state so> each mode M G TM starts out 
in the state (initM,SM), where initM is the non-default entry point of M and 
Sq[M.V] = sm- Note that as long as there is a mode M whose control state is at 
initM, no continuous steps are possible. However, any discrete step of such mode 
will come from and bring the control state of M to dx. Therefore, any 

execution of an agent A = (TM,V,I) with \TM\ = k will start with exactly 
k discrete initialization steps. At that point, every top-level mode of A will be 
at its default exit point, allowing an alternation of continuous steps from R'^ 
and discrete steps from R(jf The choice of a continuous step involving all 
modes or a discrete step in one of the modes is left to the environment. Before 
each discrete step, there is an environment step, which takes the control point 
of the chosen mode from dx to de and leaves all the private variables of all 
top-level modes intact. After that, a discrete step of the chosen mode happens, 
bringing control back to dx. Thus, an execution of A with \TM\ = kis & sequence 

So — ^ . . . S}z — ^ . . . such that 

— for every 0 < i < k, there is M G TM such that {si[M.V], Si+i[M.y]) G 
That is, the first k steps initialize the top-level modes of A. 
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— for every i > k, one of the following holds: 

— such that / is defined on [0,t] and f{t) = Si+i, and for ev- 
ery mode M G TM, {si[M.V], f[M.V]) G M.R^; that is, the step is a 
continuous step, in which every mode takes part; 

— SjAsi+i such that for every mode M G TM, Si[M.Vp] = Si+i[M.Vp]; 
that is, the step is an environment step; 

— with i > k, there is M G TM such that {si[M.V], Si+i[M.V]) G 
M.R^^ da;i that is, the step is a discrete step by one of the modes. 

Note that environment steps in agents and in modes are different. In an agent, 
an environment step may contain only discrete steps, since all agents participate 
in every continuous step. The environment of a mode can engage in a number 
of continuous steps while the mode is inactive. 

Definition 5. (Trace semantics for agents) A trace of A is an execution of A, 
projected onto the set of its global variables. The denotational semantics of an 
agent consists of its set of global variables and its set of traces. 

Let A be a primitive agent and {init, so)-^{dx, Si)^(c2, 52)^ . . . s„) 

be a trace of its top-level mode. It is easy to see that so-^si^S2^ • • • is a 

trace of A. A similar statement is true for agents with multiple top-level modes. 

4.3 Operations on Agents 

Variable hiding. The hiding operator makes a set of agent variables private. 
Given an agent A = (TM,V,I), the agent A\{ 14 } = {TM,V' , I) with V( = 
ViUVh,V( = Vg — Vh. A trace of A, projected onto the set of global variables of 
A\{V/i}, is a trace of A\{ 14 }. 

Variable renaming. Variable renaming replaces a set of variables in an agent 
A with another set of variables. Let V\ = {a;i,... ,x„},V2 = {vu--- ,J/n} be 
indexed sets of variables with Vi C A.V. Then, A[Vi := V2] is an agent with 
the set of global variables {A.Vg — V\) U V2- Semantics of the variable renaming 
operator is given by renaming the variables in the traces of the agent. 

Parallel composition. The composition of the two agents A1HA2 is an agent 
A = (TM,VJ) defined as follows:A.TM = A^.TM U A2.TM, A.Vg = Ai.Vg U 
A2.VgjA.V1 = Ai.Vi U A 2 .V 1 , and if s G A.I then s[Ai.V] G Ai.I and .s[A 2 .V] G 
A 2 .I . 

5 Compositionality Results 

We show that our semantics is compositional for both modes and agents. First, 
the set of traces of a mode can be computed from the definition of the mode 
itself and the semantics of its submodes. Second, the set of traces of a composite 
agent can be computed from the semantics of its sub-agents. For the lack of 
space, we omit the proofs and concentrate on intuitions for the results. 
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Fig. 3. Compositionality rules for modes 



5.1 Compositionality of Modes 

In order to show that our trace semantics for modes is compositional, we need 
to be able to define the semantics of a mode only in terms of the semantics of 
its submodes. 

Compositional Trace Construction. First, we show that every trace of a 
mode can be constructed using the traces of the submodes. 

Theorem 1. The set of traces of a mode M can be computed from the set of 
traces of its submodes, its closed transition relation c(T) and the set of con- 
straints Cons. 

Theoremd relies on the following observation. Given a submode N of M, we 
can “project” a trace a of M onto N and obtain a trace of N. This projection 
will 1) restrict all data states and flows to the global variables of N, 2) replace 
every subsequence of a where N is inactive into a single environment step, and 3) 
convert continuous steps of M into continuous steps of N by removing transitions 
from N.dx to M.dx and from M.de to N.de. The critical point in proving this 
observation is that, whenever the control state is at dx of M, and N is the 
active submode of M, N has its control state at N.dx, since only default exit 
transitions and the identity transition of the mode can end at dx. 

Mode Refinement. The trace semantics leads to a natural notion of refinement 
between modes: a mode M refines N if it has the same global variables and 
control points, and every trace of M is a trace of N. 

Definition 6. (Refinement) A mode M and a mode N are said to be compatible 
if M.Vg — N.Vg, M.E=N.E and M.X=N.X . Given two compatible modes M 
and N, M refines N, denoted if LMf=Lj^ . 

For a finite index set J, we write {Mi | i G /} ^ {Ni | i G /} if Mi < Ni 
for each i G L The refinement operator is compositional with respect to the 
encapsulation: 

Theorem 2. (Submode compositionality) Given a mode N, suppose SM < SN 
and let M = N[SM/ SN]. Then M N. 

The refinement rule is explained visually in Figure 0 left. If we consider a sub- 
mode N within a mode M , the remaining submodes of M and the transitions 
of M can be viewed as an environment or mode context for N. In other words. 
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a context for Ni . . . Nk is a mode M[Gi, . . . Gk] with holes or most general sub- 
modes Gi, l<i<k that have the same interface as Ni, have no local variables and 
put no constraints on the update of global variables. Two contexts are said to 
be compatible if they are compatible as modes and they also are compatible on 
their holes. 

Definition 7. (Context traces) An execution of a mode context G with holes 
Gi . . . Gk is a path 

(eo) so)^(ei) si)— I . . . — ^(e„, s„) 

through the graph of TZ of G with Xi = e for each Ci, e^+i such that Ci is in G.X 
and Ci+i is in G.E or et is in Gj.E and Cj+i is in Gj.X, for l<j<fc. A trace 
of G is obtained by projecting an execution on its global variables. 

As with modes, the set of traces of a context C is denoted by Lc and refinement 
is defined by language inclusion. Given a context G with holes Gi, . . .Gk and 
a set of modes Ni, . . . Nk such that Ni ^ Gi for l<i<k, we write G[A^i, . . . Nk] 
the mode obtained by filling the holes Gi of G with Ni. Contexts are also com- 
positional. 

Theorem 3. (Context compositionality) Let C\ and C 2 be compatible contexts 
with holes Gi . . . Gfc. If C\ :< C 2 then Gi[A^i, . . . , Nk] ^ G 2 [iVi, . . . , Nk] for any 
set Ni,l<i<k of modes compatible with the holes, i.e., Ni ^ Gi for all i. 

A visual representation of this rule is shown in Figure 0 right. The compo- 
sitionality rules allow us to decompose the proof obligation into refinement of 
submodes in the most general context, and refinement of contexts under the 
most general submode. 



Normal 




Normal ’ 


( „ , ) 




( ) 


ff e X \ 

)t^0 

t=W \dx de jl 




r 'V» 

t<10 \dx de y 


C * ° ) 




C * ° ) 


d{t) = 1 

(level s [2,10]j 




d(t) = 1 

(level G [2,10]} 



Fig. 4. Refinement example 



Consider mode Normal in Figure [D as a two-place context. Let Normal’ 
differ from Normal only by allowing rate computation to happen more often. 
The transition to Compute has a relaxed guard t < 10, as shown in Figure E| 
By Theorem 0 Normal [Maintain, Compute] A Normal’ [Maintain, Compute]. If 
Controller ’ is the agent in which Normal ’ replaces Normal, then by TheoremEI 
Controller A Controller’. 
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5.2 Compositionality of Agents 

An agent is, in essence, a set of top level modes that interleave their discrete 
transitions and synchronize their flows, the compositionality results for modes 
lift in a natural way to agents too. The operations on agents are compositional 
with respect to refinement. 

Definition 8. (Refinement) An agent A and an agent B are said to be com- 
patible if A. Vg = B.Vg. Agent A refines a eompatible agent B, denoted A<B, if 
LaQLb- 



Theorem 4. (Agent compositionality) Given compatible agents such that 
A<B,Ai<Bi and A2<B2- Let Vi = {a;i,... ,a;„},V'2 = {yi,--- ,2/n} be in- 
dexed sets of variables with Vi C A.V and let Vh Q A.V. Then A\{t 4 } ^ 
B\{Vh},A[Vi := ^2] ^ B[Vi := and A1WA2 < B1WB2 

In our example, TEink| jController ^ Tank||Controller’ by Theorem^ 

6 Conclusions 

We have presented a hierarchical modular semantics for hybrid systems. The 
proposed semantics is compositional both with respect to the system architec- 
ture (parallel agents and their subagents) and the system behavior (modes and 
their submodes). We have introduced the notion of refinement between the sys- 
tem components - both modes and agents - and showed that, in the proposed 
semantics, composition of components preserves refinement. 

We are currently working to build upon the presented compositionality re- 
sults and provide assume-guarantee proof rules for hybrid systems, extending 
the results of Pj. The proposed semantics have been used in the modeling lan- 
guage Charon | 3 | and its toolkit, currently under development by the authors. 
For further details, see 

http : / /www . cis . upenn . edu/mobies/charon/. 
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Abstract. We consider an optimal-reachability problem for a timed au- 
tomaton with respect to a linear cost fnnction which results in a weighted 
timed automaton. Onr solution to this optimization problem consists of 
reducing it to a (parametric) shortest-path problem for a finite directed 
graph. The directed graph we construct is a refinement of the region au- 
tomaton due to Alur and Dill. We present an exponential time algorithm 
to solve the shortest-path problem for weighted timed automata starting 
from a single state, and a doubly-exponential time algorithm to solve 
this problem starting from a zone of the state space. 



1 Introduction 

Timed automata !XD94| are widely accepted as a formalism to model the be- 
haviour of real-time systems: a discrete transition graph is equipped with a finite 
set of clock variables which are used to express timing constraints. Automated 
analysis of timed automata relies on the construction of a finite quotient of the 
infinite space of clock valuations. In particular, this construction is suitable to 
perform reachability analysis. Given two states s and f of a timed automaton A, 
the reachability problem can be stated as the problem of determining if there 
exists a run of A from s to t. Reachability is a core problem in system verification 
and directly applies to the verification of safety properties. 

In the theory of timed automata there are many decision problems which are 
undecidable, and decidability is in general hard. In this paper we are interested in 
an optimal-reachability problem for timed automata. Time-optimal reachability 
was first considered in , where the problem of computing lower and upper 

bounds on time delays in timed automata was solved. Minimal-time reachability 
is also considered in fNTYOOI . In |AGHh,S| . a weight w is associated with each 
location q such that w gives the cost of a unit of time spent in q. Then, given 
a cost interval I and two states s and t, the decision problem “is t reachable 
from s at a cost c G 17” {duration-bounded reachability) is addressed and solved. 

* This work is partially supported by the DARPA/ITO MoBIES grant F33615-00-C- 
1707, the NSF Career award CCR97-34115, the SRC award 99-TJ-688, the MURST 
grant TOSCA, the DARPA JFACC grant N66001-99-C-8510, and the University of 
Pennsylvania Research Foundation. 
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Here we solve a more general optimal-reachability problem, that has been inde- 
pendently solved also in | |BHF~*~| . We consider weighted timed automata, that is 
timed automata with weights (different costs) on both locations and transitions. 
The cost of a run is given by the sum of costs of the taken switches plus the sum 
of the costs associated with the visited locations multiplied for the time spent in 
each of them. Our optimization problem, which we call optimal-run problem, can 
be formalized as a tuple containing a weighted timed automaton, a source zone 
and a target zone. If the source zone contains only a state of the automaton, we 
refer to this problem as the single-source optimal-run problem. 

Our solution to the optimal-run problem consists of two main steps: first we 
reduce the optimal-run problem to a shortest-path problem in directed graphs, 
then we solve the latter. The first step is obtained by constructing a finite graph 
which is a refinement of the region automaton imm . Each clock region is split 
into several disjoint subregions relatively to a starting state and to sequences 
of resets that may occur in “potential” optimal runs. This construction is pa- 
rameterized on the differences of two consecutive fractional parts from the clock 
valuation of the starting state. When we consider a general source zone, we leave 
unspecified these parameters and the above construction reduces the optimal-run 
problem for weighted timed automata to a parametric shortest-path problem in 
directed graphs. We give a fix-point computation algorithm to solve this prob- 
lem, so obtaining a doubly-exponential time algorithm solving the optimal-run 
problem. In case the input automaton has only one clock variable, this result can 
be improved to a single exponential by adapting to our case the algorithm given 
in fK()S1fYT()9l j for solving a particular case of parametric shortest-path prob- 
lem. In case the source zone is a singleton we substitutes the parameters with 
the actual values from the starting state, and thus our optimization problem 
is reduced to a standard shortest-path problem. Using Dijkstra’s algorithm, we 
obtain an exponential time algorithm for the single-source optimal-run problem. 

The optimal-reachability problem is strictly related to other decision prob- 
lems, and in particular to the problem of synthesizing an optimal controller. 
The optimal- control synthesis problem can be informally stated as the prob- 
lem of designing a control which is able to drive, at a minimum cost, the sys- 
tem into a given target zone. In the literature, control synthesis problems have 
been considered in the context of discrete automata lUhuh'ilThohRI . timed au- 
tomata linear hybrid automata [WW7] . and general 

hybrid systems jLTS99ISPSQQj . The design of an optimal control for hybrid sys- 
tems is not trivial and in general is undecidable. The approach presented in 
this paper, can be adapted to solve the optimal-control synthesis problem for 
weighted timed automata. We observe that this generalizes the results obtained 
in !XM99| on the synthesis of a time-optimal controller for a timed automaton. 

The rest of the paper is organized as follows. In section El we define the 
optimal-run problems and we give some examples. In section El we introduce 
a graph construction to reduce the optimal-run problems to the corresponding 
shortest-path problems in directed graphs. In section 0 we present our solutions 
to the single-source optimal-run problem and to the general case. 
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2 Preliminaries 

In this section we define the single-source and the parametric optimal-run prob- 
lems. We start introducing some notation and the definition of timed automaton. 

Given a set C of n variables, a fc-zone is a subset of K" that can be obtained 
as a boolean combination of inequalities of the form x<y + c, x<y + c, x<c, 
and X < c where x,y G C and c G {0, 1, ... , k}. We denote by True the clock 
constraint which is true for any clock values. We denote by Z{C) the set of all 
the A:-zones, for all k G N. A function A : K" — >■ R." is called a reset function if 
it is equal to the identity on some of the coordinates and zero on the others. We 
denote by the set of all reset functions over K". A timed automatoi^ A is a 
tuple (Q, G, Z\, Inv) where: 

— Q is a finite set of locations; 

— G is a finite set of n clock variables; 

— Z\ is a finite subset of Q x Z{C) x A„ x Q; 

— Inv : Q — > Z{C) maps each location q to its invariant Inv(q). 

A state is a tuple (g, v) where q G Q and u G M”. We denote by S' = QxR" the 
set of states for A. A discrete step is (g, v)-^{q' , v') where e = (g, A, g') G A, 

V satisfies (5, v' = A(u), and v' satisfies Inv(g'). A time step is (g, u) — ^ 
where v' = v + t, t > 0, and ly + t' satisfies Inv(g) for all 0 < t' < t. 

A step is (g, u)-^(g', i/') where (g, u) — ^ (g, u") and {q,v”)-^{q' ,v'), for 
some v” G R", that is a transition e taken after spending some time t in 
the current location. A run r of a timed automaton A is a finite sequence 

(<?o,uo)4^ (gi,ui)-^ . . (gfc_i,Ufe_i)^ (gfc,Ufc). We say that r starts at 
(901 uo) and ends at (g*,, u^). The definition of r allows time to be spent after tak- 
ing the last transition Ck-i- A weighted timed automaton is a timed automaton 
A with the following cost functions: 

— Js ■ A — > N {switch cost), and 

— Jd '■ Q — > N {duration cost). 

Given a run r of A and cost functions Jg, and Jd, associate costs to r as 
follows: 

— Js{r) = J2i=i Js{ei), and 

“ Jd{f) = X)i=0 ■ Jd{gi)- 

The total cost associated to a run r is then J{r) = J s{r) + Jd{f)- We are inter- 
ested in determining optimal-cost runs for a timed automaton. In the following 
examples we informally introduce some notions that we will formalize in the rest 
of the section. 

^ The standard definition of timed automata requires also an acceptance condition 
and a symbol alphabet. Since we are not interested in studying languages accepted 
by timed automata we omit these features here. 
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Example 1. Consider the timed automaton defined in Figure^such that Jd(0) = 
3, Jd(l) = 1) and the switch costs are all 1. Suppose that we start from state 
s = (0,x,y) for 0 < x,y < 2 and we want to reach a state in location 2. 
Possible minimal-cost runs from s to a state s' = {2,x',y') are either ri = 

{0,x,y)~^ {l,xi,yi)^ {2,x + 2-y,2), or X 2 = {0,x,y)-^ {2,2,y + 2-x) 
for ^3 = (2 — x) (obviously, staying in location 2 longer might only increase the 
overall cost). According to the cost function J, the cost of ri is Js{ri) + Jd{ri) = 
2+3ti+{2—y—ti) = 4—y+2ti and the cost of V 2 is Js{r 2 )+Jd{r 2 ) = H-3(2— a;) = 
7 — 3x. Clearly, J{ri) is minimized when t\ = 0, that is the transition from 0 to 1 
is taken immediately. Moreover, assuming ti = 0, J(ri) < J(r 2 ) if y > 3(x — 1), 
and J{ti) > J{r 2 ), otherwise. Thus, a minimal-cost run from s to a state in 
location 2 depends on the clock valuation of state s. 




Fig. 1. A timed automaton with more than an optimal run from a same location. 



Example 2. Consider the timed automaton defined in FigureElsuch that Jd(0) = 

1, Jd(l) = 2, and the switch costs are all 1. Suppose that we start from 
state s = (0, x) for 0 < X < 2 and we want to reach a state in location 

2. Possible minimal-cost runs from s to a state s' = (2,x') are given by 

rt = (0,x)-^ (l,xi)-^ (2,2). Notice that vt is a run parameterized by t, 
where t is the time at which the first edge is taken. Thus Jir^) = J s{rt)EJ d(ji) = 
2-|-t-l-2(2 — t — x) = 6 — t — 2x. Hence the cost of is minimized if t is maximized. 
Since t < {2 — x) must hold, the optimal cost for a run starting at s is (4 — x), 
but none of the runs starting at s has such a cost. In fact, for any actual run rt 
there exists a ^ > 0 such that t= {2 — x — S,), and J(xt) = (4 — x-h^). Vice-versa, 
for any ^ > 0 there exists a run r such that J(r) = (4 — x -I- ^). Clearly, there 
is not a minimal-cost run but we can determine a run whose cost is arbitrarily 
close to the optimal one. 

Now we formalize the notion of optimal cost, optimal run, and approximation 
of an optimal run. Given a timed automaton A, a state s, and a target zone T, 
an optimal cost for a run from s to T is a J* such that J* < J(r) for any run r 
from s to a state in T, and for any ^ > 0 there is a run r such that J(r) < J* -1-.^. 
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Fig. 2. A timed automaton with no optimal runs from a location. 



If there exists a run r* such that J(r*) = J*, then r* is said to be an optimal 
run. As shown in Example 0 sometimes an optimal run from a state s to a 
target zone T does not exist. In these cases, we are interested in a family R 
of runs such that all the runs coincide on the sequence of switches and for any 
^ G M_i_ there exists a run r G R such that J{r) < J* +^, where J* is the optimal 
cost over all runs from s to T. That is we can determine a sequence of runs in 
R whose costs are arbitrarily close to J* . We call such a family of runs R an 
approximation of an optimal run. Given a timed automaton A, a source zone S, 
and a target zone T, we consider the problem of determining an optimal run from 
a given state s G S' to T, if one exists, or an approximation of an optimal run, 
otherwise. We call this problem a single-source optimal-run problem. We also 
consider a more general problem, a zone optimal-run problem, defined as the 
problem of determining a symbolic representation of the solution to the single- 
source optimal-run problem for all states in S. In Example Q if we consider as 
target region all the states in location 2 and as only source state (0,0,0), then a 
solution to the corresponding instance of the single-source optimal-run problem 
is ri with = 0. As observed in Example ^ if we consider as source zone the set 
of states {0,x,y) such that 0 < a;,y < 1, then the solution of the corresponding 
instance of the zone optimal-run problem is r\ with t\ — 0 ii y > 3(a; — 1), and 
C 2 , otherwise. 

We end this section with an example on an air-trafhc control problem that 
we will use subsequently in the paper. 



Example 3. Consider the timed automaton in Figure 0 It models a scenario in 
which two aircraft send a landing request to an airport, and our goal is to allow 
both the aircraft to land safely and at minimum cost. Safety requires that only 
one aircraft at a time must be acknowledged for landing, thus there are two 
possible choices: aircraft 1 waits for the landing of aircraft 2 to be completed, 
or vice-versa. There are costs c\ and C 2 to pay for forcing respectively aircraft 
1 and aircraft 2 to wait. Moreover, there is also a cost, expressed by Wi, which 
is related to the time spent waiting. Alternatively, aircraft i can make, at a cost 
c(, a maneuver that allows to spend w[ instead of wt per each time unit. This 
maneuver takes at least time 1. Since it is realistic to reduce the time a runway 
stays unused, we penalize this event by a cost cq per time unit. Finally, we 
assume that the landing of each aircraft takes at least time 1 since the related 
acknowledgement was issued by the control tower. 
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Fig. 3. An air-traffic control problem. 



3 The Graph Construction 

In this section we give the graph construction underlying the reduction of the 
single-source optimal-run problem to the shortest-path problem and the zone 
optimal-run problem to a parametric shortest-path problem. The obtained graph 
is a refinement of the region automaton of a timed automaton, in the 

sense that each vertex v carries more information than a region. This additional 
information mainly concerns the sequence of resets needed to reach v from a 
starting vertex, and the construction preserves the transitions of the region au- 
tomaton. Via this construction we emphasize the states of the timed automaton 
that might be visited in some optimal runs. We start by recalling the concepts 
of labelled directed graph and region automaton, then we describe our graph 
construction. 

Let 6> be a set of real-valued parameters, we denote by D the set of linear 
expressions over 0. Given an alphabet S, a D-lahelled directed graph G is a pair 
(V,E), where V is a set of vertices, and E C V x D xV is a set of 15-labelled edges. 

A path 7T from vq to Vn in G is a sequence vq v\ u„_i Vn 

f 1 

such that Vi-i Vi G E for i = 1, . . . ,n. For a path tt, the cost of tt is given by 
Sr=i /*■ ^ path TT from v to v' is a shortest path if tt is the path with minimum 
cost among those connecting v to v'. Notice that varying the values of parameters 
in 0 the shortest path of a graph may change, that is to different valuations of 
parameters may correspond different sets of shortest paths in the graph. 

Consider now a timed automaton A. By definition its set of states is infinite. 
However, they can be partitioned in a finite number of equivalence classes, called 
regions, which are defined by a location and a clock region. Denoted by Cx the 
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largest constant in clock constraints involving the clock variable x, a clock region 
is described by: 

— a constraint of type c — 1 < x < c, x > c^, or x = c for each clock variable 
X and c < Cx', 

— the ordering of the fractional parts of the clock variables x such that x < Cx- 

Thus a clock region denotes a set of clock valuations. Given a clock valuation 
V, [v] denotes the clock region containing v. A state {q, v) belongs to a region 
(q',a) if q = q' and v € a. A clock region a is said to be open if for any clock 
variable x and c < Cx, x = c does not hold in a. Otherwise a is said to be 
a boundary clock region. These definitions apply to regions in an obvious way. 
The key property of this equivalence, is that all the valuations belonging to a 
region satisfy the same set of clock constraints from the given timed automaton. 
Consistently we say that a clock region a satisfies a constraint <5 if satisfies S 
for any n G a. A clock region a' is said to be a time-successor of a clock region 
a if and only if for any v G a there is a, d G such that v -\- d G a' . The region 
automaton of A is a transition system defined by: 

— the set of states R{S) = {{q, a) \ q G Q and a is a clock region for A}; 

— the transition rules R{A) such that: {{q,a) , {q' , a')) G R{A) if and only if 
{q, A, 6, q') G A and there is a time-successor a" of a such that a" satisfies 5 
and a' = [A — >■ 0]a". 

We denote the region automaton corresponding to A as i?(A). For the sake of 
simplicity, in the following when no confusion can arise we refer to the value of 
a clock variable x by a; itself. With x we denote the fractional part of a clock 
variable x. Let s = {q, v) be a state of A and (0 «2 ■ • ■ x'n ~iv+i 1) 

be the ordering of the fractional parts of the region containing a clock valuation 
V (notice that is either = or <). With §{s) = (di, . . . ,'diq+i) we denote the 
differences between consecutive values in the above ordering, that is "di = x^, 
^N+i = 1 ~ 2 ;(y, and di = x' — x[_i for i = 2, . . . ,N. In the following we 
will use (-di, . . . , dAT+i) to denote these differences in the starting state. The 
graph we are going to define is parameterized over (-di, . . . , dAr+i). Moreover, 
for i,j < N, we denote by I{i,j) the set of integers {i, . . . ,j — 1}, if t < j, and 
{i,. ■ ■ ,N} U {1, . . . , j — 1}, otherwise. 

The region automaton does not carry enough information to solve our op- 
timization problems. Thus we define a labelled directed graph whose vertices 
correspond to “sub-states” of the region automaton. For a given state {q,a') 
of the region automaton, a sub-state {q, a) is such that a is a convex region 
contained in a' . Denoted by (0 x'^ «2 • • • x'h 1) the ordering of the 

fractional parts in a clock region a' , we consider sub-regions a of a' such that 
for some of the «i’s which are equal to <, the difference between xi_i and x( is 
very close to 0. Thus we represent a by a' and specifying in the ordering of the 
fractional parts if a < is relative to a “small” difference (denoted by or to a 
“large” difference (denoted by <). We call each such sub-region a a boundary 
sub-region. Intuitively, the reason we are interested in boundary sub-regions is 
that the cost functions we consider are linear, and their infimum over a given 
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region is reached on the boundary. Thus optimal runs leave open regions from 
states which are arbitrarily close to their boundaries. As a consequence optimal 
runs visit also states characterized by having clocks values either with arbitrarily 
close fractional parts or with fractional parts which reflects the starting state and 
the reset history of the computation. For this reason, we add to each boundary 
sub-region a tuple of indices (A, . . . , ik) from {1, . . . , n -I- 1} such that: k is the 
number of large differences in the ordering of the fractional parts, i/ corresponds 
to the /-th large difference in the ordering of the fractional parts, and there exists 
a, d G {1, . . . , A:} such that id+h < id+h+i for h = 0, . . . ,k — 1, where the sums 
{d + h+ 1) and (d -I- h) are modulo k. We call such tuples distance tuples, since 
they are used to store the difference between two consecutive fractional parts 
when this difference is “large” (i.e., they are not arbitrarily close). We define the 
set of vertices V as the set of tuples {q, a, (ii, . . . , ik)) where g is a location, a is 
a boundary sub-region, and (A, . . . , ik) is a distance tuple from {1, . . . , n -I- 1}. 
For a vertex (q, a, (A, ... , ik)}, the sum A) gives the time to leave the 

region since this subregion is entered. 

The set of edges B contains three types of edges: immediate switches, time 
edges and delayed switches. Informally, immediate switches correspond to tran- 
sitions taken in the current state, time edges correspond to letting time elapse 
until the next region is reached, and delayed switches correspond to transitions 
taken at the “beginning” or at the “end” of the closest open region (this region 
if it is an open region, the next otherwise). 

Given two vertices v = {q, a, (zi, . . . , ih)) and v' = {q', (3, (ji, . . . ,jk))), there 
J s{e) 

is an immediate switch v v' if there exists a transition e of R{A) from {q, a') 
to {q' , j3'), where a' and P' are respectively the regions of R{A) containing a and 
P, and the sequence (ji, . . . ,jk) is obtained from (ii, . . . , ih) by deleting all the 
indices ii such that all the clocks between the /-th and the (/ -I- l)-th large 
differences (in the ordering of the fractional parts of a') are reset in e. 

Consider a vertex v = {q,a,{i\, . . . , ih)) and let (0 |/i «2 • ■ • Vk 

1) be the ordering of the fractional parts in a. If we assume that ot{yk) -I- 1 is not 
larger than the largest constant in the timing constraints involving yk (i.e., when 
time elapses the first integer value reached by yk is at most this constant), we 
add to if a time edge v — ^ v' for v' = {q, P, {ji, . . . ,jh')) where P is the closest 
time-successor of a such that the conditions expressed by one of the rows of the 
following Table 1 are satisfied (where (0 y'l «2 • ■ • ~k v'k ~k+i 1) denotes 
the ordering of the fractional parts in P, and I = 2, ... ,k): 
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In the other case, time edges are defined in the same way except for the 
fact that the clock yk does not appear in the ordering of the fractional parts of 
v' since it has reached its highest constant. To see an example of a time edge. 
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consider a vertex v = {q, 0<x<y<z<l, (1, 2, 3, 4)). By row 1 of the above 
table we have a time edge from v to {q, 0<x<y<lAz = l, (4, 2, 3)). The 
distance tuple (4,2,3) captures the fact that time {1 — z) has elapsed and thus 
the distance in time from a; to 0 is increased by (1 — z), the fractional part of z 
is now 0, and all the other distances stay unchanged. 

Given a vertex v £ V as above, we add to if a delayed switch v — ^ v" for 

J s{e) 

any vertex v" £ V such that there exists an immediate switch v' v" and 
c = c' + Js(e), where v' = {q, /3, (ji, . . . ,jh')) and (3 is the closest time-successor 
of a such that the conditions expressed by one of the rows of the following Table 
2 are satisfied (where (0 y'^ ■ ■ ■ ~kV'k ^k+i 1) denotes the ordering of the 

fractional parts in /3, and I = 2, . . . , k): 
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For a given tuple of parameters = (di, . . . , we denote by Ga(i?) the 

ii-labelled directed graph {V, E). We recall that for our purposes li represents the 
differences between the fractional parts of two consecutive clocks in the ordering 
of the fractional parts in the starting state. The construction of Ga{'&) is general 
in the sense that it does not depend on the particular source and target zones of 
the problem, but only on the timed automaton. This allows us to use it for solving 
both the single-source optimal-run problem (for a fixed d) and the zone optimal- 
run problem (d belongs to a convex set). As an example of application of the 
above construction, we discuss a fragment of the graph Ga{3)) for the weighted 
timed automaton modelling the air-traffic control problem from Example 0 (see 
Figure 0). For the sake of simplicity, we have marked with 1, ... ,5 the vertices 
of Ga{3)) in Figured, and we refer to them by these numbers. Consider vertex 
1. Since in the timed automaton from Figure 0 there is a transition from W\ to 
W[ resetting clock x\, we have in Ga(^) an immediate switch from 1 to 2. Edges 
from 1 to 3 and from 1 to 4 are delayed switches obtained by the same transition 
above and respectively rows 3 and 4 of Table 2. The edge from 1 to 5 is a time 
edge and is defined by row 2 of Table 1. Notice that for a given state s = (q, v), we 
have corresponding vertices of Ga{’3{s)) of form (q, a, (A, • ■ • , *fc)), where v £ a. 
Moreover, each edge is labelled by the actual cost of the corresponding “activity” 
in A, that is for immediate switches we have just the cost of the A transition, 
for time edges the cost of spending the time upto the end of the current region 
in the current A location, and for delayed switches the cost corresponding to the 
A transition plus the cost for the time spent in the current location before that 
the transition is taken. We have the following lemma. 



Lemma 1. Given a timed automaton A, the size o/G^('d) is exponential in the 
length of clock constraints of A. 
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Fig. 4. A fragment of Ga(i?) for the weighted timed automaton in Example 0 



Proof. In [A I )H4j the authors proved that the size of the region automaton is 
0(|A| where |i5(A)| denotes the length of the clock constraints. A sim- 

ple counting argument gives that the number of ways to substitute < with < 
in the ordering of the fractional parts of a clock region is at most and 

the number of tuples of indices we use to represent the relative differences be- 
tween the fractional parts is at most n2". Thus the size of Ga{'&) is at most 
0(|A| n2^"+^ and since n = 0(|5(A)|), it is exponential in the length of 

the clock constraints. 

4 Optimal-Runs in Weighted Timed Automata 

4.1 Single- Source Case 

In this section we prove that the single-source optimal-run problem in timed 
automata can be reduced to the shortest path problem in a weighted directed 
graph. To see this we introduce first some notation. Let sq be a state {qo,vo) 
of a weighted timed automaton A and 'd(so) = denote by 

g{so) the vertex (go, UO) (*o,i) • ■ • j *o.Afo)) Ca(i^(so)) such that vq G oq and 
ioj is the j-th largest distance in the ordering of the fractional parts in oq- 
Given a positive real f « 1 and a path tt = (^o? (*o,ii ■ ■ ■ i *o,Afo)) 

{qt, ah, {ih,i, ■ ■ ■ ,ih,Nh)) in Ga(^?(so)), we 
denote by i? 7 r(C) tbe set of runs of A starting at sq and obtained by replac- 
ing with {qj,Vj)~^ {qk,i'k) each portion {qj,aj,{ijA, ■ ■ ■ ,ij,Nj)) ^ 

{qk,ak, (ifc.i, • • ■ Gfe.TvJ) of tt such that: 

- {qj-i,aj-i,{ij-iA, ■ ■ ■ ^ {qj,aj,{ijA, ■ ■ ■ ,ij,Nj)) is either an 

immediate or a delayed switch; 

-for I = j...,k - 2, {qi,ai,{iiA,...,ii,Ni)) ^ 

(q/+i,a;+i, (b+ 1 , 1 , . . .,ii+i^Ni+J) is a time edge; 
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- • ■ ■ ,ik-i,Nk-i)) {qk,ak, (*fc,i> • ■ • is either an 

immediate or a delayed switch. Let tj = t' + r" and vj + t' € ak-i- In the 
case of an immediate switch t" = 0, while in the other case r" is such that: 

— if the delayed switch is obtained by rows 1 and 2 of Table 2, then Vj+tj G 
Ofe-i and the largest fractional part in Vj + tj is greater than (1 — ^); 

— otherwise, denoted as a' the time-successor of ak-i which is first entered 
by letting time elapse from a valuation in ak-i, it holds that Vj+tj G a' , 
moreover if the delayed switch is obtained by rows 4 and 5 of Table 2, 
the largest meaningful fractional part in Vj + tj is greater than (1 — ^), 
and if the delayed switch is obtained by rows 3 and 5 of Table 2, the 
smallest meaningful fractional part in Vj + tj is less than 

— ej is the transition corresponding to {qk-i,ak-i, • ■ • ,'ik-i,Nk-i)) 

{qk, 0!k, {ik,l, • ■ • , * fc , Affc ))- 

In the following we assume that ^ is a positive real number such that ^ << 1. 
By the definition of Ga{'&) and i? 7 r(C), we have the following lemma. 

Lemma 2. Given a timed automaton A and a state s = {q, v) of A, if tt is a 
path of GA{'d{s)) from g{s) of cost Ct^ then R-k{£.) is a set of runs of A such that 
for any e > 0 there exists an r & Ruif) such that Ct^ < J{r) < Ct^ + e. 

To complete our reduction we need the following lemma. 

Lemma 3. Given a run r of A from a state s to a target zone T, there exists 
a path TT of GA('d{s)) from g{s) to a vertex corresponding to a state in T such 
that the cost of tt is not larger than J(r). 

Proof. The interesting case is when transitions in r are from states that do not 
belong to any of the subregions encoded by Ga(i^(s)) vertices. Assume that A 
in run r takes a transition e from an open region a after spending some time in 
it, and e is the first transition in r with this property. Clearly, upto e, r has a 
corresponding path tt in G^(i9(s)) whose cost is not more than J{r). We observe 
that by definition there must be two delayed transitions e\ and of Ga{'&{s)) 
corresponding respectively to the cases e is taken as soon as a is entered and 
e is taken just before leaving a. Moreover, consider two A runs r\ and V 2 that 
differ from r only for the fact that in ri A takes e after an arbitrarily short time 
spent in a, while in r 2 A takes e after an arbitrarily short time before leaving 
a. Clearly, J(r) > min{ J(ri), J(r 2 )} holds. Thus we can add to tt the transition 
corresponding to the run with the least cost between ri and r 2 - Applying 
iteratively this argument, we determine a path tt in GA(d(s)) of cost c < J{r). 

As a direct consequence of Lemmas 0and0 we have the following theorems. 



Theorem 1. Given a timed automaton A, a state s of A, a target zone T, tt 
is a shortest path of GAi'&{s)) starting from g{s) to a vertex corresponding to a 
state in T if and only if is an approximation of an optimal run of A from 

s to T. 
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Theorem 2. Given a timed automaton A, a state s of A, a target zone T, there 
exists an optimal run of A from s to T if and only if for a shortest path tt of 
G^(r?(s)) from g{s) to a vertex corresponding to a state in T there exists a run 
r G such that the cost of tt is equal to J{r). Moreover, r is an optimal 

run of A from s to T . 

Given a timed automaton A, a source state s, and a target zone T, the 
following algorithm solves the single-source optimal-run problem: 

1. Let G be the graph obtained from GA('d(s)) by collapsing all the vertices 
corresponding to a state in T in a single vertex vt- 

2. Solve the single-source shortest-path problem on G from g{s). 

3. Let 7T be a shortest path from Vg to Vt- OutpulH and the cost of tt. 



Theorem 3. The single-source optimal-run problem can be solved in time ex- 
ponential in the size of the timed automaton. 



4.2 The Algorithm for the General Case 

In this section we consider the zone optimal-run problem. We give an exponential 
time algorithm to solve this problem for timed automata with at most 1 clock 
and a fix-point algorithm in doubly-exponential time, for the general case. 

We start considering the general case. Since we want to solve the problem 
of determining the optimal runs from any state of the source zone S' to a state 
of a target zone T, for parameters d in we consider only values given 

by -d = d{s) for a state in s G S. Thus it holds that -di -I- . . . -I- = 1 

and we can eliminate a parameter by the substitution = 1 — 

From now on, we will assume that d{s) is the tuple (i?i, . . . and G^(t?(s)) 
is the graph obtained after the substitution 'dAr+i = 1 — The algo- 

rithm that we are giving, labels the vertices of G^('d) with sets of linear ex- 
pressions on d = {di, . . . The meaning of these expressions is that given 

a state s G S the minimum over these expressions gives the optimal cost of 
a run from s. An expression is a first-degree polynomial in d\, . . . ,dN , and 
(1 — with integer coefficients. That is, an expression has the form 

f{d) = 00 - 1 - aidi -I- ... -I- aNdN + OAr+i(l - ^0: where oq, . . . , oat+i are 

nonnegative integer constants. We denote expressions by {N -\- 2)-tuples of co- 
efficients and write (oq, . . . , oat+i) for the above expression f{d). We denote by 
^ the natural extension to tuples of the total ordering < over reals. Moreover, 
let /, /' be two expressions, and v,v' be two vertices of GA(d), (f,v) -< {f ,v') 
if and only if / ^ /'. A set A of tuples of type (/, u), for an expression / and a 
vertex v, is said to be minimized (with respect to if for any (/, v), (/', v') G X, 
(f,v) and if',v') are not comparable with respect to A. 

^ This step needs a further refinement to distinguish between an approximate solution 
and an optimal solution. It is not entirely straightforward, but it can be handled at 
the same complexity. We defer the reader to the full version of the paper. 
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The algorithm we present computes a labelling function I that maps any 
vertex u of C?y 4 ('(^) to a minimized set of pairs (/, v) for which there exist a path 
7T and a state s € S such that: 

— 7T is a shortest path of Ga{^{s)) from it to a vertex corresponding to T, 

— the first edge e of tt connects u to v, and 

— the cost of TT is given by /(i?(s)). 

We can summarize our algorithm in the following steps: 

1. Initialize I by assigning l{u) = {(0, . . . , 0, u)} for u corresponding to a state 
in T, and /(it) = 0 for all remaining vertices. 

2. repeat 

I' ^ 1; I ^ UPDATE(r) 

until I' = I 

3. Output 1. 

We just need to specify the function Update. Consider an edge e a vertex 
It = , ifi))- We have the following cases: 

— e is an immediate switch from it to v: for (ag, . . . , UAr+i, v') G define 

(og, . . . , u) such that Og = oq + Ce, and a' = for i = 1, . . . , (N + 1), 
where Ce is the cost of e; 

— e is a time edge from it to v: for any (og, . . . , ajv+i, u') G l'(v), define 

(flg, . . . , such that if e is obtained by rows 1 and 2 of Table 1 and 

1 G then a' = + Jd{q), otherwise a' = af, 

— the edge e is a delayed switch from it to v: for any (og, . . . , atv+i, 'cO ^ 
define (a'g, . . . , a(y_|_^, u) such that if e is obtained by rows 1, 2 and 4 of Table 

2 and i G then o' = + Jd{q), otherwise o' = o^. 

Let /"(it) be the set of all the tuples generated for it. After executing / ^ 
Update(/'), /(it) contains the set obtained deleting from /'(it) U /"(it) all the 
tuples (/, u) such that f f for some (/',u') G l'(u) U l"(u). Moreover, once 
the function / is output, it is easy to determine the optimal cost and generate 
the corresponding solution from / and the graph G^(d), given d. We observe 
that each of the tuples (/, u) belonging to l{u) corresponds to a path from u to 
a target vertex. Thus the cardinality of /(it) is bounded above by the number of 
simple paths in Ga(i?). Hence we have the following theorem. 

Theorem 4. The zone optimal-run problem can be solved in doubly- exponential 
time. 

If we restrict to timed automata with just one clock variable, it is possible to 
solve the zone optimal-run problem in singly exponential time. We consider the 
algorithm given in |K()81IYT()91j to solve a particular shortest-path problem 
with only a parameter and edge costs given by (c — i?), for constants c. This 
algorithm runs in polynomial time and can be modified in order to obtain a 
polynomial time algorithm to solve the parametric shortest-path problem with 
edge costs given by a first-degree polynomial of r? (z? S [0, 1]). 

Theorem 5. The zone optimal-run problem for automata with one clock vari- 
able can be solved in exponential time. 
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Abstract. Reach set computations are of fundamental importance in 
control theory. We consider the reach set problem for open-loop systems 
described by parametric inhomogeneous linear differential systems and 
use real quantifier elimination methods to get exact and approximate 
solutions. The method employs a reduction of the forward and back- 
ward reach set and control parameter set problems to the transcendental 
implicitization problems for the components of special solutions of sim- 
pler non-parametric systems. For simple elementary functions we give an 
exact calculation of the cases where exact semialgebraic transcendental 
implicitization is possible. For the negative cases we provide approximate 
alternating using discrete point checking or safe estimations of reach sets 
and control parameter sets. Examples are computed using the redlog 
and QEPCAD packages. 



1 Introduction 



Today integrated systems which combine physical processes with information 
systems {i.e. digital programs) are in great demand. In fact complex systems 
which have been designed recently incorporate both differential equations to 
model the continuous behavior and discrete event systems to model instanta- 
neous state changes in response to events. Systems that are finite state machines 
with differential equations at each discrete state are called Hybrid Systems. 

A lot of research effort has been devoted to develop mathematical models, 
specification formalisms, analysis/design/control methods and tools to help con- 
trol engineers in building such systems (see \ 1 81,41 )l2b] ) . Most of the applications 
of hybrid systems are safety critical. Safety is usually encoded as avoidance of an 
undesirable region of the state space. Consequently, the most important prob- 
lems for analyzing hybrid systems are verification problems; these are essentially 
reachability problems, that ask whether trajectories of the hybrid systems reach 
certain undesirable (unsafe) regions from an initial region. 

Computing the reach set of hybrid systems is difficult because hybrid systems 
have an infinite state space. Due to the difficulty of computing the reach set for 
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systems of differential equations, formal verification methods and tools for hybrid 
systems have been developed fM7\ . These methods and tools, however, can deal 
with only very simple continuous models as, e.g. x = 1, Ax = b. What is actually 
required is to handle hybrid systems with more complicated continuous parts. 

Decidability of reachability problem for hybrid systems with linear differen- 
tial equation of the form y = Ay + Bu is discussed in ISE3. This is a significant 
class of linear differential equations that is widely used in linear control theory. 
The results are based on the notion of “o-minimality” m from model theory 
and “quantifier elimination” El. O-minimality is used to define a class of hy- 
brid systems “o-minimal hybrid systems” and it is shown that all o-minimal 
hybrid systems admit finite bisimulations in |22ll . To make the bisimulation al- 
gorithm computationally feasible, they utilize mathematical logic, in particular, 
real quantifier elimination, as main tool to represent and manipulate sets sym- 
bolically. Since quantifier elimination, in general, is possible for the polynomial 
theory of reals they have found subclasses of o-minimal hybrid systems that 
are definable in the theory. 

Remark: There are many results that apply quantifier elimination to control theory 
ymm . In IL!Si;-il quantifier elimination is used for verification problems (reachability 
and observability problems) of discrete-time polynomial systems. 

In this paper we study in particular reach set problems for continuous open- 
loop systems described by parametric systems of linear differential equations El. 
Roughly speaking reach set problems are concerned with the relations between 
possible values of the state variables at some initial time to the corresponding 

values at later points in time. The specific problems studied in this paper are 
the following: 

1. Fix a set M of values of the state variables at to; what are the possible 
corresponding values at later points t in time (up to some bound ti or oo). 
(Forward reach set) 

2. Fix a set N of “safe” values of the state variables. Find a set M as large 
as possible of initial values of the state variables at time to that guarantees 
that the values of the state variables will for all later time points t (up to 
some bound t\ or oo) remain inside N. (Backward reach set) 

3. Fix a set M of values of the state variables at to and a set N of “safe” 
values of the state variables. Find a set P as large as possible of the control 
parameters such that all state variables with initial values at to in M will 
have values in N for all later time points t (up to some bound ti or oo.) 
( Control parameter set) 

Our main tool is the method of real quantifier elimination in computer algebra. 
This approach was introduced into reach set computations in m- In a series of 
papers they showed how to get exact solutions of the forward reach set problem 
for certain homogeneous linear differential systems of special type with constant 
coefficients El and for associated inhomogeneous systems with very special 
right hand side m- The exact solutions are always obtained as semialgebraic 
sets described by a boolean combination of polynomial inequalities. 



Reach Set Computations Using Real Quantifier Elimination 



65 



Here we extend this ad hoc approach for special types of differential systems 
to a systematic study of the type of results obtainable by an approach via real 
quantifier elimination. By reducing the approach to its bare essentials, we obtain 
a much wider systematic framework applicable to a considerably larger class 
of systems. The main observation is that all the problems mentioned above 
can be reduced by exact symbolic algorithms to an implicitization problem for 
certain basic transcendental functions associated with the given system. Exact 
solutions for implicitization problems with rational parametrizations are well- 
known . Here we deal with the corresponding problem for transcendental 
parametrizations that has been studied only for special cases e.g. in H3I2DI. 

Our main results are as follows: We associate with every parametric linear 
system of differential equations y = A{t)y + b{t,r) a finite system F of basic 
functions. Then for semialgebraic sets M, N all three problems can be solved 
exactly by real quantifier elimination relative to the implicitization problem for 
the components of the functions in F. Moreover the discrete point version of 
these problems require only finitely many evaluations of functions in F. We 
prove a theorem that determines the exact classes of vector- valued functions of 
the kind arising in linear differential systems with constant coefficients, where 
exact semialgebraic implicitization is possible. As a corollary we obtain the exact 
limitations of the approach of for linear differential systems with constant 

coefficients and special right hand sides. 

We propose several ways to overcome these limitations by approximate com- 
putations: One way is to compute exact reach sets at a finite selection of discrete 
time points. This is always possible and practically quite efficient, but may lead 
to underestimation of the true forward reach set, depending on the selection of 
time points. Another approach separates the common time variable into differ- 
ent time variables. This leads to an overestimation in the implicitization problem 
resulting in an overestimation of the forward reach set and an underestimation 
of the backward reach set and the control parameter set: So all three approxi- 
mations are on the safe side. 

We illustrate some problems and solution methods by examples computed 
in the REDLOG-package of reduce M and QEPCAD m- We expect that our 
results can be extended to the hybrid systems with linear continuous parts. 



2 Reach Sets and Transcendental Implicitization Problem 

2.1 Problem Statement 

We consider parametric inhomogeneous systems S of linear differential systems 
of the form y = A(t)y + b(t, r) with an n x n matrix A(t) of real continuous 
functions aij(t) and a vector- valued real continuous function b(t,r) defined on 
some interval /. The inhomogeneous part is assumed to be a linear combination 
b{t,r) = X^i=i with continuous functions gi : I — > K", and real parame- 

ters Ti . Such a system can be viewed as an continuous open-loop control system 
with control parameters r = (ri, . . . ,rfc). Let M be some subset of R." and fix 
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an initial time point to € I : Then we denote the set of all solution functions 
/ : I — > R" of the given system with parameters r= (ri, . . . ,rk), by TV, and 
the set of all solution functions f G with initial value /(to) G M by F = Fm,t- 
We consider the following forward reach set problems: 

discrete reach sets Compute for finitely many time points t\ <...< tm ba I 
the union of the sets {f{ti) \ f G Fm^t}- 
bounded reach set Compute for a given time ti > to in I the set {/(t) | / G 
FM,nto ^ t ^ ti}- 

unbounded reach set Suppose I 3 [to,oo), and compute the set {/(t) | / G 
Fm,t, to < t}- 

All computations should be performed in explicit dependence on the control 
parameters r. Any solution of the discrete reach sets problem yields an lower 
estimate for the sets to be computed in the bounded and unbounded reach set 
problems. 

Of equal interest are the corresponding “backward” reach set problems that 
are a kind of “dual” to the corresponding “forward” problems. 

Some backward reach set problems are as follows: Let be a subset of R". 

backward discrete reach sets Compute for finitely many time points ti < 
. . . < in / the sets {/(to) I /(ti), . . . , /(tm) e N}. 
backward bounded reach set Compute for a given time t\ > to in J the set 
{/(^o) I /(t) e N for all to < t < h}. 

backward unbounded reach set Suppose / O [to,oo), and compute the set 
{/(^o) I /(t) e N for all to < t}. 

From the viewpoint of control theory these problems have still other vari- 
ants concerning the determination of suitable control parameter values r = 
(ri, . . . , Tfc). Let M as before be a subset of R”, and let N be another subset of 
R". Then we have the following natural control parameter set problems: 

discrete point control Compute for finitely many time points ti < . . . < tm 
in I the set {r G R^ | f{ti) G N for all / G Fm,v, 1 < * < m}- 
bounded interval control Compute for a given time ti > to in / the set 
{r G R'^ I /(t) G N for all / G Fm,t, to < t < ti}. 
unbounded interval control Suppose / 3 [to,oo), and compute the set (r G 
R^ I /(t) G N for all / G ^ ^}- 

In order to make these problems mathematically precise, we need to specify 
the way in which the input sets M and N, and the output sets should be de- 
scribed. For an approach using symbolic computations it is natural to consider 
semialgebraic sets as possible inputs. These are subsets of R" described by a 
boolean combination :p{xi, . . . , Xn) of real polynomial inequalities. If in addition 
all the polynomials involved in :p{xi, . . . ,Xn) are linear, then the set described 
by ip is called semilinear 

Our goal is to solve the forward and backward reach set and control parame- 
ter set problems for semialgebraic input sets as far as possible with descriptions 
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of semialgebraic sets as outputs. This, however, is not always possible. Hence we 
consider also the computation of overestimating the forward reach sets and un- 
derestimating the backward reach set and the control parameter sets by suitable 
semialgebraic sets. 

Our main tool will be a reduction of reach set and control parameter set 
computations to corresponding implicitization problems for a fixed finite sys- 
tem of functions associated with S, namely a fundamental system /i , . . . , /„ for 
the homogeneous system Sq associated with S, and special solutions hi of the 
parameter-free inhomogeneous system Si given by y = A{t)y+gi{t) for 1 < i < k. 
We refer to {/i, . . . , /„, hi , . . . , hk} as a system of basic functions for S. 

Implicitization problems for rational parametrizations of algebraic varieties 
have been widely considered in computer algebra PP7] . Here we have to study 
the corresponding problem for the vector- valued functions /i, hi, . . . ,hk, 

arising from the system S. As these functions will in general be transcendental, 
we refer to these problems as transcendental implicitization problems. 

More precisely, we consider the following transcendental implicitization prob- 
lems for given functions fi : I — > K" for 1 < * < fc : 

discrete points implicitization Compute for finitely many time points t\ < 
. . . < tm CO- I the values (/i(ti), . . . fk{ti)), regarded as points in 
bounded implicitization Compute for a given time ti > Iq in I the set 
{{flit ), . . . fk{t)) e K”'' \to<t< ti}. 
unbounded implicitization Suppose I A [to,oo), and compute the set 

{(/i(t), ...^(^))eR"'=|^o<n• 

The first problem amounts to simple evaluations of the given functions. No- 
tice that the unbounded and bounded implicitization problem for a single solu- 
tion of the differential system S is in fact a special case of the unbounded and 
bounded forward reach set problem for S, respectively, namely for the case of a 
singleton set M. 

2.2 Reduction to Implicitization Problems 

Next we show that all reach set computations and control parameter set com- 
putations listed above can for semialgebraic input sets M, N be reduced in an 
exact symbolic way to one of these implicitization problems. All these reductions 
require real quantifier elimination as fundamental tool. For the case of discrete 
points forward and backward reach set and control parameter set and semilinear 
input sets M, N we find moreover that the output sets are also semilinear. 

Let ip{xi, . . . , Xn) and ifixi, . . . , Xn) be quantifier-free formulas describing 
the semialgebraic input sets M and N, respectively. Let y = Ay -|- h{t, r) with 
t>{t,r) = ’^rigiit) be a parametric linear system S with control parameter 
Let fi be a fundamental system of solutions of y = Ay. Let hi be a special 
solution of the system y = Ay + gi{t). Then by the superposition principle, a 
special solution of the system S is given by X/i=i Note that here rfs may 
be regarded as constants or as free parameters. Then it is straightforward to 
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write down first-order formulas describing the respective forward and backward 
reach sets and control parameter sets in terms of evaluations of the basic func- 
tions /i, . . . , /„, hi, , hj^, the given formulas i.p{xi, . . . , a;„), ip{xi, . . . , Xn) and 
a quantifier-free formula ^(yn, . . . , yin, ■ ■ ■ , Vni, ■ ■ ■ , Vnn) describing the com- 
bined range of (/i, . . . , fn, hi, . . . , hk), as a semialgebraic set. All these formulas 
will involve several quantifiers over real numbers. By real quantifier elimination 
one can construct equivalent quantifier-free formulas, and thus get the desired 
semialgebraic descriptions. 

We will exhibit concrete first-order formulas for some reach set problems and 
control parameter set problem. The remaining cases are handled similarly in 
p]. The forward discrete reach set problem can be described by the following 
formula and hence be solved by real quantifier elimination and evaluation of the 
basic functions at finitely many points. 

3x1 . . .3xnMJ2i^ifi + T,ir^h^)ito) A [A"=i % = (Ei + E* 

V ■ • • V /\j-i Vj = (Ei ^ifij + Ei fihij){tm)]). 

Next suppose we have a quantifier-free formula y{yii, . . . , yin, ■ ■ ■ , yni, ■ ■ ■ , 
ynn, Zii, ..., Zin, Zki, ..., Zkn) describing the combined range of (/i , . . . , fn, 
hi, . . . , hk) on the interval [to, oo) or [to, ti[. So /i(yn, . . . , Zkn) holds for n{k+n)- 
tuple in if and only if this tuple is in the combined range of (/i, . . . , 

hi, . . . , hk) on the given interval. Then the forward bounded and unbounded 
reach set problem, respectively, can be described by the following formula and 
hence solved by real quantifier elimination: 

3xi . . .3Xn[ip{J2i^ifi + J2i^i^i)ito) A 3yu . ..3ynn3zii . ..3Zkn{h{yil,- ■ -,Zkn) 

A Aj = l ttt = + Ei'^’i^d))]- 

With the same formula y,, the backward bounded and unbounded reach set 
problem, respectively, can be described by the following formula and hence solved 
by real quantifier elimination: 

3xi . . . 3x„[A”=i Vj = (Ei + Ei nhj)(to)) A Vyn . . . Vy nn V.2:ii . . . ^ Zkn 
{y{yii,...,Zkn) -t 1|^{T,^x^yi + T,^n^){t))]■ 

Finally, the bounded interval control problem and the unbounded interval 
control problem, respectively, can be described by the following formula and 
hence solved by real quantifier elimination: 

3xi . . .3Xn[^{J2i^ifi + Ei^i^i)(^o) A Vj/n . . .VynnVzii . ..Vzkn 
{y{yii,...,Zkn) -t A(E*2^i^+ E*^i3i)(^))]- 

As a corollary to these semialgebraic parametric descriptions of reach sets 
we also obtain semialgebraic descriptions of the corresponding reach sets, where 
the control parameters range over a prescribed semialgebraic set C. 

Corollary 1. Let C C K.* be a semialgebraic sets of control parameters de- 
scribed by a quantifier-free formula 7 (ri, . . . , r^). Let p{yi, . . . , yn, xi,. . . , Xk) be 
a quantifier-free formula describing a forward /backward reach set relative to the 
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control parameters r. Then the corresponding reach set for arbitrary control pa- 
rameter values in C is described by the formula 3 ri . . . 3rfc(7(r) A p{y,r)), and 
Vti . . . Vrfc(7(r) — >■ p{y,r)), respectively, and hence is also a semialgebraic set. 



Example 1 . Consider the inhomogeneous system Q y = Ay + b with 



A := 



1 -3 

-2 2 



b := ri 



,2t 



r2 



„4t 



0 



Then the basic functions are 2 ® J . A quantifier- 

free formula p{yn, yi2, 2/21, ?/22, -^ii, -^12, -^21, Z22) describing the combined range of 
these functions for t G [ 0 , 00) is obtained as follows: Notice that the range of 
on [ 0 ,oo) is exactly ( 0 , 1 ], and that e* = = l/(e“‘)^,e^‘ = l/(e“‘)^. 

So /i can be taken as the formula 

< 1 A 3yi2 = 2yii A 2/212/11 = -1 A 2/222/11 = 1 A 

- n A — 1 A — 



0 < 2/11 < 1 A 3 yi 2 = 2 yii A 2/212/11 = -1 A 2/222 
zii = 0 A 3 zi 22 /n = 1 A 2 z 2 iyii = l ^ Z22 



= 0 . 



3 Exact Transcendental Implicitization 

Here we consider cases, where the unbounded and bounded transcendental im- 
plicitization problem for given functions fi : I — >■ K." (1 < i < A:) has an 
exact solution. Notice that the transcendental implicitization problem refers 
only to the component functions fij(t) of fi(t); the grouping of these compo- 
nent functions into vector- valued functions is irrelevant here. So we may as- 
sume w.l.o.g. that k = 1 and that we deal with a single vector-valued function 
f(t) := (/i(t), . . . , fn{t)). Then the exact transcendental implicitization problem 
is to determine the range of f(f) on an unbounded interval [to, 00), or a compact 
interval [to,fi] contained in I. Since the / is continuous, this range is always a 
connected subset of M" . 

In particular for n = 1 the range is a real interval J; moreover J is compact 
for the bounded implicitization case. In the unbounded implicitization case J is 
compact iff / is bounded on [to, 00), otherwise it is a closed semiinfinite interval 
or all of K.. In particular J is always a semialgebraic set that can computed 
explicitly from upper and lower bounds for /. In other words the unbounded 
and the bounded transcendental implicitization problem always has a positive 
solution for n = 1 . 

For n = 2 there are two well-known cases, where exact unbounded and 
bounded implicitization is possible, namely the sin-cos-pair and the sinh-cosh- 
pair: If / has components /i := cos{p{x)),f2 ■= sin(p(x)), where p{x) is a real 
polynomial of positive degree, then the range of p{x) on [to, 00) includes an 
unbounded interval; consequently the range of / on [to, 00) is exactly the unit 
circle {(a;i,a;2) \ xf x^ = 1 }. On a bounded interval [to,ti], the range of p{x) 



^ This is taken from |1()| (p.586 example 3.13). 
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is again a compact interval, and so the range of / is a connected subset of the 
circle that can be easily computed as semialgebraic set from the range of p{x). 
For the hyperbolic case, where /i := cosh(p(cc )), /2 := sinh(p(a:)), the situa- 
tion is analogous, except that the role of the circle is replaces by the hyperbola 
{{xi,X 2 ) \x\-xl = 1}. 

The next theorem shows that exact transcendental implicitization is pre- 
served under composition of functions in a very general sense: 

Theorem 1. Let f{t) := (/i(t), . . . , fk{t)) he a vector valued function such that 
the range of f on every compact or unbounded closed interval I is a semialge- 
braic set described by a quantifier-free formula ipj{xi, . . . ,Xk)- Let g be a con- 
tinuous real function defined on some compact or upper semiinfinite closed in- 
terval L' . Let hi (1 < i < n) be semialgebraic real functions defined on some 
subset o/ K" extending the range of f. Let pi{xi, . . . ,Xn,y) be quantifier-free 
formulas defining the graph {{xi, . . . ,Xn,y) \ y = hi{xi, . . . ,Xn)} of hi. Then 
the vector-valued function f*{t) := (/*(t), . . . , /*(t)) with components f*{t) := 
hi{fi{g(t )), . . . , fn{g{t))) for 1 < i < n has a semialgebraic range described by the 
formula 'tjj{xi,...,Xn) ■= 3yi . . . 3j/„((/?j(2/i, . . . , y„) A Pi{yi, ■ ■ ■ ,yn,Xi)), 

where J is the range of g(t) on L'. 

The proof is obvious. Notice that the algorithmic quantifier elimination for the 
ordered field of real numbers this formula is required in order to transform the 
formula if into an equivalent quantifier-free formula that describes the range of 
/* as a semialgebraic set. Typical instances of g and hi are real polynomials or 
real rational functions. The method can in particular be applied to the situation, 
where / consists of a sin-cos-pair or a sinh-cosh-pair as described above. Other 
interesting examples are pairs (p, p'), where p(<) is a Weierstrass p- function 
IP. Then the range of (p, p'), on a large enough interval is a real elliptic curve 
{(x, y) \ y^ = 4x^ — g 2 X — gs}. See jS| for the more examples. 

4 Semialgebraic Implicitization for Simple Elementary 
Functions 

In this section we characterize those cases of linear differential systems S with 
constant coefficients and “simple right hand side”, where an exact implicitiza- 
tion of the system of basic functions for S is possible. The condition on the right 
hand side b{t) of the system is as follows: All components bi{t) of b{t) are R-linear 
combinations of functions of the form t‘^* cos{u>it) , sin(wit), where di 

are non-negative integers and ai,uji, o.i are real numbers. Then it is well known 
that a special solution of the inhomogeneous system and the fundamental solu- 
tions of the homogeneous system are again real linear combinations of functions 
of this kind. We call linear systems of this form regular and functions of type 
cos(wt), sin(o;f), with a,uj,a real numbers simple elementary func- 
tions. In some special cases of regular systems, it has been shown how to solve 
the reach set problem by an implicit semialgebraic implicitization of functions 
of the following type in : (i) real polynomials Pi{t), (ii) exponential 
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functions with rational values of ai. (iii) trigonometric functions cos{u>it), 
sin(o;it), for rational u>i. 

In the following we show that for simple elementary functions there are only 
few more cases which allow unbounded exact semialgebraic implicitization; all 
these cases are covered by Theorem ^ of the last section. In most of the re- 
maining cases the exact semialgebraic implicitization problem is unsolvable. In 
fact we provide a complete characterization of those cases, where unbounded 
semialgebraic implicitization is possible. 

Let f(t) := (/i(t), . . . , fn{t)) with non-constant, pairwise different component 
functions fi{t) := cos(wit), or fi{t) := sin(wit), where di are non- 

negative integers and ai,tUi are real numbers. Moreover we assume that the 
functions fi appear in cos-sin-pairs, whenever yf 0. 

Theorem 2. Let f : — > K" be as above and let n > 2. Then the range 

of f is a semialgebraic set ijf one of the following holds: 

1. For all 1 < i < n, fi(t) := . 

2. For all 1 < i < n, di = 0, fi{t) := and dimQ{span{ai , . . . , a„)) < 1. 

3. For all 1 < i < n, di ^ 0, 0, fi{t) := and 

dimQ{span{ai , . . . , a„)) < 1, and ^ 

4- For all 1 < i < n, fi{t) := cos(wit), or fi{t) := sin(o;it), and 
dinvQ{span{uJi, . . . ,Wn)) < 1- 

Moreover in these positive eases a quantifier-free formula describing the range 
of f can be computed algorithmically over the reals. 

Idea of the Proof. In the cases mentioned above the unbounded semialgebraic 
implicitization is always achieved by the methods of the previous section, in 
particular Theorem E It remains to show that in all other cases the range of / 
is not a semialgebraic set. This requires a case distinction. In each case we show 
that the assumption that the range of / is semialgebraic leads to a contradiction. 
Based on the assumption that the range of / is semialgebraic we construct new 
semialgebraic sets with impossible properties. Either this set is one dimensional 
such that neither the set nor its complement is a finite union of intervals or it 
describes the graph of a semialgebraic function with an impossible rate of growth 
(compare P|). See [3| for details of the proof. 

This theorem clearly shows the limitations of the approach presented in m 
ED. In fact we have the following immediate corollary: 

Corollary 2. Let y = Ay with constant n x n-matrix A be a homogeneous 
system of linear differential equations. Then exact semialgebraic implicitization 
is possible for a fundamental system of solutions of the system iff one of the 
following cases holds: 

1. All eigenvalues of A are zero, i.e. A is a nilpotent matrix. 

2. All eigenvalues Ai,...,A„ of A are non-zero, pairwise distinct reals, and 
dimQ{span{Xi , . . . , A„)) < 1. 

3. All eigenvalues Ai, . . . , A„ of A are purely imaginary, say of the form Ai = 

with non-zero pairwise distinet reals yn, and dimQ{span{iii , . . . , /in)) < 1. 
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5 Approximate Solutions 

In this section we study the cases, where an exact semialgebraic unbounded or 
bounded implicitization is definitely not possible. In these cases we want to find 
a semialgebraic superset of the true forward reach set and a semialgebraic subset 
of the true backward reach set or the true control parameter set, both if possible 
such that the set difference to the true reach set or control parameter set is in 
some sense “small enough.” Then an inspection of the reduction formulas shows 
that an overestimation of the implicitization problem leads to an overestimation 
of the forward reach set and an underestimation of the backward reach set and of 
the control parameter set z.e. for “safe” estimations. Hence we are reduced to the 
problem of finding a semialgebraic superset of the true range of a transcendental 
vector valued function on a compact or upper semiinfinite closed interval. 

One strategy to find overestimations of the range is separation of variables: 
It comes in two flavours: Separation of variables in different components, and 
separation of variables in products. 

Separation of variables in different components : Let f{t) = (/i(t) 

,...,/„(t)) be defined on an interval I. Then separation of variables in dif- 
ferent components yields the function g{t) = , fn{tn)) defined on the 

cube /" with range((/) O range(/). The range of g is easily computed as a box 
Ji X • • • X J„, where Ji is the range of fi. Notice that this box is in fact the 
smallest box containing the range of /. 

Separation of variables in products : Suppose the component functions 
of the given functions are products fi{t) := /i,i(t) • • • fi,m{t), where each 
is defined on the interval I. Put gj{t) := {fij, ■ ■ ■ , fn,j)'^ ■ Then each gj is also 
defined on the interval I. Let Bj be the range gj, and put C := 
where the multiplication is performed on the elements componentwise. Then C 
is obviously a superset of the range of /. 

Example 2. Let / be the upper semiinfinite interval [0, oo). 

1. Let /i := cos(t), f^ '■= sin(t). Then the true range of / is the unit circle. 
Separation of variables in different components yields as overestimation the 
closed unit square. 

2. Let /i := cosh(t), fi := sinh(t). Then the true range of / is the hyperbola 

{{x,y) \ = 1}. Separation of variables in different components yields 

as overestimation the “quadrant” {{x,y) \ x,y > 1}. 

3. Let /i := e‘cos(t), /2 := e*sin(t). Then the true range of / is an expanding 
exponential spiral. Separation of variables in different components yields as 
overestimation the full plane Separation of variables in products yields 
as better overestimation the annulus {{x,y) \ x'^ + y'^ > 1}. 

4. Let /i := cos(t), /2 := e~* sin(t). Then the true range of / is a contracting 
exponential spiral. Separation of variables in different components yields as 
overestimation a closed box {{x,y) \ —e^ < x < 1, — < y < 
Separation of variables in products yields as overestimation the closed disk 
{{x,y) I x“^ + y^ < 1}. These approximations are incomparable. So their 
intersection is a common improvement of both. 
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6 Complexity 

In this section we briefly discuss the complexity of our algorithms. From the 
results on complexity of quantifier elimination in Pj we can give upper bounds 
for the asymptotic complexity of our approach: 

Discrete point reach set problems are described by purely existential formu- 
las. Hence the complexity of quantifier elimination is at most simply exponential 
in the dimension of the differential system. For fixed dimension it the compu- 
tation runs in a polynomial time. The complexity of bounded and unbounded 
reach set problems is the same as for the discrete reach set problem for a fixed 
number m of points. The backward discrete reach set problems can be solved in 
singly exponential time. The complexity of backward bounded and unbounded 
reach set computation is of type ' (generalized singly exponential). The up- 
per complexity bounds for the control parameter set problems are same as for 
the corresponding backward reach set problems. 



7 Computational Example in redlog and qepcad 



In this section we report on experimental results in reach set and control param- 
eter set computation. In 0 we have presented experimental results for numerous 
examples that illustrate the different problem types and solution methods. Here 
we display only one of these examples with non-constant coefficients to show 
the generality of the approach. All computations are performed in the redlog 
package m of REDUCE 3.7 and qepcad mE . The main algorithm employed is 
the linear and quadratic quantifier elimination P^TTT) of REDLOG and quantifier 
elimination based on cylindrical algebraic decomposition m of QEPCAD. 



Example 3. Consider the inhomogeneous system y = Ay + b with 



A ■= 





b := ri 



/ 2t cos(t^) \ 
2f sin(f^) J 



Then basic functions are ^ system we 

illustrate the computations in the forward/backward unbounded reach set and 
the control parameter set problems below (Note that we set to = 0): 

• Forward unbounded reach set: A quantifier-free formula y{yu,yi2, ?/2i, 2/22, 
-2^11, -2^12) is obtained from the following first-order formula y,o 

fXo = 3 u 3 v{u^+v'^ = IA2/11 = vAyi2 = uAy2i = uAy22 = —vAzu = vAz\2 = 0 )) 



by using quantifier elimination. By using REDLOG we have 

M := 2/11 -I- 2/?2 - 1 = 0 2/11 + 2/22 = 0 A 2/11 - zii = 0 A 2/12 - 2/21 = 0 A Z12 = 0 

^ All the computations are executed on a SUN SPARC station Ultra I (140MHz). 



74 



H. Anai and V. Weispfenning 



in 10 ms. Then we set ri = 1 and moreover </^ = (0 < a:i < 1 A X 2 = 0)- 
Then forward unbounded reach set problem is solved by using real quantifier 
elimination for the following first-order formula /reach; 

/reach = 3x\{ip A /reachaux) 

where 

/reachaux = 3yll3yl23y2l^y22^zu_3zl2{^J. /\yi= xiyn + X 22/21 + 

A?/2 = Xiyi 2 + X 2 y 22 + ^’1^12) 

By using QEPCAD for /reach we obtain as an answer for the forward unbounded 
reach set; yf + 4y| — 4 <= 0 in 10 ms. 

• Backward unbounded reach set: /r is the same formula as in forward 
unbounded reach set. We also set ri = 1 and tp{xi,X 2 ) = (—5 < Xi < ^ A — | < 
Xi < |). Then the backward unbounded reach set problem is solved by using 
real quantifier elimination for the following first-order formula breach; 

breach = 3xi3x2{yi = X 2 A y 2 = x\ A breachaux) 

where 

breachaux = VyiiV?/i 2 V?/ 2 iV 2 / 22 V 2 :iiVzi 2 (/r -)► (-| < Xiyn + X 22/21 + ^ 

A - 5 < Xiyi 2 + X2J/22 + rizi 2 < 5) 

By using REDLOG for breach we obtain in 420 ms a semialgebraic description of 
the backward unbounded reach set consisting of 21 atomic formulas. 

• Control parameter set: The formula y, is the same as in the reach set cases. 
We also set (^ = (0 < Xi < 1 A a ;2 = 0) and ’4i{x\,X2) = {—^<xi<^A—^< 
xi < ^). Then control parameter set problem is solved by using real quantifier 
elimination for the following first-order formula pcontrol; 

control = 3xi{ip A controlaux) 

where 

controlaux = V?/iiVj/i 2 Vy 2 iVy 22 VziiVzi 2 (/r -t (-| < Xiyn + X 22/21 + riZn < | 
A - i < Xiyi2 + X2y22 + riZi2 < 1) 

By using REDLOG for control we obtain in 70 ms a semialgebraic description of 
control parameter set consisting of 12 atomic formulas. It can be simplified to 
the result — 1 < ri < | by hand calculation. 

8 Conclusions 

In this paper we have studied forward and backward reach set and control pa- 
rameter set problems for continuous parametric open-loop systems described by 
a system of parametric linear differential equations with arbitrary coefficients. 

The approach using quantifier elimination was introduced into reach set com- 
putations in m- We extend their ad hoc approach for special types of differential 
systems to a systematic study of the type of results obtainable by an approach 
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via real quantifier elimination. Thus we obtain a much wider systematic frame- 
work applicable to a considerably larger class of systems. The main observation 
is that all the problems can be reduced by exact symbolic algorithms to an im- 
plicitization problem for certain basic transcendental functions associated with 
the given system. 

We have proved a theorem that determines the exact classes of vector- valued 
functions of the kind arising in linear differential systems with constant coeffi- 
cients, where exact semialgebraic implicitization is possible. As a corollary we 
have obtained the exact limitations of the approach of for linear differ- 

ential systems with constant coefficients and simple elementary inhomogeneous 
part. We have also proposed several ways to overcome these limitations by ap- 
proximate computations. The problems have been illustrated by examples com- 
puted in the REDLOG-package of reduce and qepcad. 

Further research will be concerned with an extension of these results to hybrid 
systems. 
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Abstract. In this work we present a novel control design methodology 
for under-actuated mechanical systems. As part of the design process we 
use the reachability analysis tool d/dt [ARI)M90BDnn| to see whether 
there is a switching sequence which can drive the system to a desired 
periodic orbit. Much of the work in the design of the control law is done 
manually using classical control te chniques (unlike the fully-automatic 
approach advocated in |ARD+nn| l. and d/dt is used to complement 
these techniques. We hope this work will contribute to the proliferation 
of reachability-based techniques to the control engineer’s tool box. 



1 Introduction 

The algorithmic approach to the analysis of hybrid systems, first put forward 
explicitly in |ACH~*~9,^ . is inspired by a computer science approach to verifica- 
tion of automata. The system under consideration is viewed as a generator of 
trajectories and the problem of verification consists of checking whether there 
is an individual trajectory which violates some specification, e.g. reaches a bad 
state. Likewise, the controller synthesis problem is phrased as restricting sys- 
tematically the set of all possible behaviors in order to satisfy a property. The 
algorithmic approach consists in making a brute-force search in the state-space, 
based only on the description of the system dynamics. Initially this approach 
has been applied to restricted classes of hybrid systems where the continuous 
dynamics has a constant derivative in every state, see e.g. jsnni] for timed 
automata, and |A(1H+9,5IAMP9,5IHHW97| for hybrid automata. More recently 
attempts have been made to lift this approach to systems with non-trivial dy- 
namics. In particular, some of the authors were involved in the development of 
d/dt, a tool for verification and controller synthesis for hybrid systems with lin- 
ear continuous dynamics |ABDM99iDQQ< . The synthesis algorithm implemented 

* This work was partially supported by the European Community Esprit-LTR Project 
26270 VHS (Verification of Hybrid systems) and the French-Israeli collaboration 
project 970maefut5 (Hybrid Models of Industrial Plants). 
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in d/dt |ABD+OOm00| suggested a very idealistic scenario for switching-based 
control: the user defines the dynamics at the various modes, as well as the con- 
trol objective, and the tool automatically generates the appropriate conditions 
for mode switching. 

This approach attempts to obtain the general-purpose flavor of discrete veri- 
fication tools and it is still very remote from control engineering practice. In the 
continuous world, every class of systems has its own special character as well as 
its corresponding mathematical tricks which are used extensively by engineers 
during the controller design process. Coordinate transformations, dimensional- 
ity reduction, simplifying assumptions or linearization cannot be captured by 
straightforward reachability analysis. 

In this paper we show how reachability-based techniques can be combined 
with more “knowledge-based” methods in order to derive control strategies for a 
non-trivial class of dynamical systems, namely under-actuated mechanical sys- 
tems. We propose a general methodology for designing controllers for such sys- 
tems and demonstrate it on a double-pendulum example. The complexity of the 
system as given initially exceeds the current capabilities of reachability-based 
tools: its dynamics is non-linear and control is done using continuous actuation. 
Moreover, the system is of dimension n while the dimensionality of the available 
control is m < n. The proposed approach to control this system by switching is 
based on the following principles. 

1 . The state-space can be transformed and partitioned via a diffeomorphism (/) 
into an m-dimensional part ci and an (n — m)-dimensional part 62- 

2 . Using standard control techniques, e\ can be controlled to zero. Given this 
control, the remaining part is a closed system which defines the dynamics of 
62 (called the Zero dynamics). 

3 . Each diffeomorphism induces a different control law for its zero dynamics 
and hence a particular “mode” for the dynamics of the the uncontrolled 
part of the system. We use a parameterized family of diffeomorphisms which 
becomes finite after discretizing the parameters. 

4 . The dynamics of 62 at each mode can be linearized around its equilibrium 
point. It is possible to choose the parameters so that the linearized system 
has periodic orbits in every mode. It should be kept in mind that the validity 
of the linear model is restricted to the neighborhood of the equilibrium. 

5 . If our goal is to reach a specific periodic orbit, we can achieve it by a sequence 
of mode switchings. At each mode, however, a different quantity is controlled 
to zero. Hence, when we switch from controlling ci to controlling e^, the 
latter should already be close to zero. This restricts the parts of the state- 
space of the 62 system where switching is allowed and leads to modeling 
the system as a hybrid automaton where the transition guards reflect these 
constraints. 

The role of d/dt is then to check whether, based on the hybrid automaton 
representation, it is possible to reach from one orbit to another by mode switching 
and how much time it takes. 
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2 Control of Under- Actuated Mechanical Systems 

2.1 Under- Actuated Mechanical Systems 

We consider the class of jointed mechanical systems without flexibilities, the 
dynamics of which is given by Lagrange equations: 

M{q)q + N{q,q) = Wr (1) 

where M is the symmetric positive definite matrix defining the kinetic en- 
ergy and N gathers generalized gravity, Coriolis and centrifugal forces; q is the 
n— dimensional vector of generalized (joint) coordinates; F includes all external 
generalized forces and W is a constant matrix. 

If we now assume that the generalized forces are only actuation torques/forces 
(i.e the system is friction- free and no other potential-based actions occur), then 
the system is called under -actuated if rank {W) < n. Without loss of generality, 

we can consider that W = i ^ ) with m < n the number of actuators. 

y ^n—mxn J 



2.2 Zero Dynamics 

Let us consider a diffeomorphism 4>- 

</'(<7)= (2) 

where C\ is m-dimensional. Then, the dynamics m projected on the constraint 
Cl = 0 is called the zero dynamics associated with cj). It is given by: 

P{q){M{q)q + N{q,q))=0 (3) 

with P = In — W{JiM~^W)~^JiM~^ the projection operator, in which 
Ji = A control objective can therefore be to bring the system to this zero 
dynamics, specified by the goal task ei = 0, and to stabilize it. Since dim (ei) = 
dim (r), all the available actuation forces/ torques have to be used for that pur- 
pose. In fact, that can be done trough partial decoupling/feedback linearization: 
it can be easily seen that using the control 

r={JiM-^W)-^{u-jiq + JiM~^N) (4) 

we obtain ei = w, assumed that JiM~^W is nonsingular. It then remains to 
specify an adequate input u which stabilizes ei, asymptotically or in finite time, 
in order to drive the system to the zero dynamics. Once reached, its motion is 
then governed by eq. 0 ), which is free, since no more control is available. In 
many cases, this free motion is a periodic orbit. The idea now is to specify such 
a periodic orbit as a final goal, recalling that we can consider the choice of (j) 
as a way to modify it. The problem addressed in the following is then to study 
the reachability of this behavior starting from given initial conditions, using a 
sequence </i, 02 ■ • ■ , he successive jumps from an orbit to another one. 
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2.3 Handling the Periodic Orbits 

Let us consider the case where m = n— 1, i.e. the zero dynamics can be expressed 
using a single coordinate denoted by x\. When the phase portrait of the system 
is a closed curve O, this periodic orbit, which characterizes the zero dynamics, 
can be uniquely specified by a pair (0, X^) where is a point on the orbit, for 
example the initial conditions. Let us assume (assumption AO) that the equation 
of O in the phase plane is of the form V{x\,xi) — 1^ = 0, the invariance being 
expressed by P = 0. 1^ is a so-called Lyapunov function. For a non-actuated 
conservative mechanical system, the natural V is the mechanical energy. Since 
it is not the case here, V can only be called by analogy the “energy” level of the 
orbit . 

Let us now consider the particular case where the set of (pi consists of func- 
tions of given analytical form depending on a fc-dimensional vector of real param- 
eters p. Then p can be considered as an auxiliary control of the system. Giving 
some bounds to the parameters and the variables, so that they range over Dp 
and Dxo , respectively, the set of all possible orbits for the system is 

0 = {0(p, X°):pG Dp X° G Dxo}. 

When V is known, the set can also be parameterized by p and V . 

The problem we address now is the following: let us define a desired behavior 
of the system as a goal orbit O*; then, given an initial orbit Oq ^ O* , can we 
reach O* by modifying p? We don’t consider here related problems of automatic 
control: existence of the orbits, active stabilization, continuous control of p, which 
will be addressed in forthcoming papers. Instead, we focus our attention on 
a discrete approach, i.e. to the questions: is there a sequence of intersecting 
orbits allowing to reach O* through jumps on the parameters and how long 
time will it take? Assuming here that these jumps are instantaneous and don’t 
disturb the overall behavior (assumption Al), we can therefore forget the effect 
of the control 0 and consider for the analysis the related set of zero dynamics 
uniquely. We are therefore led back to a problem of reachability analysis of 
a hybrid system: each discrete state is an homogeneous differential equation 
associated with given values of the parameters; transitions are allowed when 
orbits of different modes are compatible with each other, i.e. when continuous 
state variables reach some particular values. We will illustrate the approach on 
the double pendulum example. 

3 The Case of the Double Pendulum 

The considered testbed is the double pendulum depicted in Figure Q] The reader 
is referred to for details on experimental issues. Terms in eq. O write 

for this system as: 



M = 



mil mi2 

mi2 TO22 



( 5 ) 
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and: 



with: 



N = 



Ni{qi,q2,qi,q2) 

^2(91, 92,91, <72) 



Cn C12 

C21 C22 



Gi 

G2 



( 6 ) 



TOll = TOiZf + 7712(^1 + Li + 2L1I2C2) 

TO12 = m2(/| + L1I2C2) 

TO22 = TO2?| 

Cii = —7772^1^252(72 

C 12 = —m2Lil2s2{qi + ( 72 ) (7) 

C21 = 7772 ^1^252(71 
C22 = 0 

Gi = -\- T7727vi)s 1 -t“ 7772^2512) 

G2 = 97772/2512 

where si := sin{qi) , C7 := cos{qi) , sij := sin{qi + qj). We consider the case 

where only the hip is actuated. Therefore W — choose the 

diffeomorphism (j) and the control F such that 

ei = qi - aq2 - b = 0 ; 62 = 92 (8) 

where a and b are two real parameter^. Therefore the zero dynamics we have 
to consider is simply: 



/ (77722 + 0^^12)92 + (0(^21 + ^22)92 + G2 — 0 
\qi = aq2 + b 



(9) 



where it assumed that 77722 + 077712 ^ 0 (assumption A2, satisfied when — < 
a < This system can be expressed in the single coordinate (72. It is a 

second order nonlinear differential equation, for which the natural state vector 



is AT = 




( 92 - 9l 
V 92 



In order to perform reachability analysis, we have to 



linearize the system. Its equilibrium points X* 




are solutions of G2(92) = 



0, i.e, for a —1 (assumption A3): 



92 = 



b + kiT 
1 + (7 



( 10 ) 



We consider in the following only the case k = 0. The equation of the system 
linearized around the center gj is: 



X = Ax = 



0 1 
—a 0 



X 



( 11 ) 



^ Note that expression 0 specifies the desired spatial trajectory of the tip of the 
double pendulum, while the “energy” level will set the amplitude and the time 
proHle of its motion along this trajectory 
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where a = I 2 + i^Licos(j^) . For ensuring the existence of periodic orbits, the 
eigenvalues of A have to be imaginary, which implies that a has to be strictly 
positive (assumption A4). The Lyapunov function associated with the system, 
i.e the energy level of an orbit is 

V=^{axl + xl) ( 12 ) 

For the purpose of reachability analysis it is more comfortable to work with 
the same system of coordinates in every state, hence we transform the linear 
dynamics of equation (HD into an affine dynamics over y = ( 92 , 92 ): 

■' = -^» + “=(-oo)»+(ay 

Finally we have to remember that the system is submitted to physical bounds on 
the joints: qi G Introducing them in (0 leads to linear constraints 

on the parameters. 

When we switch from (j) to 4>' there might be a transient period until the 
system settles in the new zero dynamics. In order to make assumption Al (tran- 
sitions are immediate) realistic we need to make sure that e\ and e\ be already 
close to their zero. For qi this means 



e'll = 111 - a' 92 - ^'1 < £i 



(14) 
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Since qi = aq 2 + b this reduces to 

\{a- a')q 2 + {b- b')\ < €i (15) 

For q 2 we need: 

|e'i| = |(a - a')(?2| < £2 (16) 

These conditions, which form rectangles in the phase-space of the zero dynamics, 
will be used as transition guards in the hybrid automaton model. Note that these 
conditions are symmetric, i.e. they are the same, in terms of <72 and q 2 for the 
transitions from {a',b') to (a, 6). Of course, their global physical interpretation 
does depend on the source state of the transition. 

The system is modeled as a hybrid automaton with 7 states, each representing 
a pair (a, 5) of parameters (Figure |5|l. At each state the dynamics is of the form 
X = Ax + u where A and u for tj^ various states are: 

Sq Si S 2 S3 S4 S5 S0 



01 01 01 01 01 01 01 
- 0.0479 0 - 0.0878 0 - 0.1167 0 - 0.1982 0 - 0.2326 0 - 0.3143 0 - 0.3555 0 

0 0 0 0 0 0 0 
0.0011 0.0000 - 0.0012 0.0000 - 0.0039 - 0.0090 - 0.0140 



The transition guards are computed according to (1151 and (II 6 H with ci = 0.05 
and C 2 = 0.02. In addition, we restrict the tranitibns toQippen between pairs 
of “close” states, i.e. |a — a'\ < 0.15 and \b — b'\ < 0.1. 




Fig. 2. The hybrid automaton for the double pendulum. The transition guards between 
pairs of states are written as products of intervals. 



84 



E. Asarin et al. 



In order to facilitate the experimentation with d/dt we have augmented the 
input syntax to include parameters and formulae referring to them. For example, 
state So and its outgoing transition is specified as: 

state: 0; 
matrixA: 

0.0 1 . 0 , 

[-12-(aO/(l+aO))*Ll*cos(bO/(l+aO))] 0.0; 

input : type convex_vert 

0.0 [(bO/(l+aO))*(-12-(aO/(l+aO))*Ll*cos(bO/(l+aO)))] ; 
transition : 
label goOl : 

if in guard: type rectauigle 

[-(-epsl+(b0-bl))/ (aO-al)] [-(epsl+(b0-bl))/ (aO-al)] , 

[eps2/ (aO-al)] [-eps2/(a0-al)] ; 
goto 1; 

4 Results 

The problem we solve with d/dt is the following: given some initial low-energy 
orbit (more preeisely, a eonneeted set of orbits) is there a sequenee of switchings 
that brings the system to its target, a higher-energy set of orbits? This problem 
is essentially a controller synthesis problem for the eventuality specification, un- 
like the safety controller synthesis that we have treated in !abd+oo| . We are 
interested in reaching the desired orbit with the least number of mode switchings. 

We illustrate informally the synthesis procedure that we employ in order to 
derive the switching controller. Consider an initial set of orbits characterized by 
the rectangle (in the (< 72 , ( 72 ) space) P = [0.7 x 0.9] x [0.01, 0.02] at state S 3 and a 
goal orbit characterized by F = [1.05, 1.3] x [0.01, 0.02] at the same state. Starting 
from the inital set (s, P) we calculate, in a breadth- first maimer , all its successors, 
i.e. continuous successors, and then, via intersection with the guards, the discrete 
successors. We continue until at some level k of the search tree, there is one or 
more paths having a leaf (s, Q) such that Q intersects F. The search graph of the 
first iteration is shown in Figure 0 and there are two intersections with the goal 
orbit after 4 transitions, along the paths S 3 , 32 , 83 , 82 , S 3 and S 3 , S 2 , Si, S 2 , S 3 . For 
every such path we do backward reachability analysis to find the predecessors 
of the goal orbit at every node and, in particular, the subset of P from which 
the goal can be reached by taking the k transitions that correspond to the path. 
This information is also used to derive the controller by restricting the guards. 
In our example we conclude that points satisfying (72 G [0.7552, 0.9] can reach 
the goal orbit by following the sequence S 3 , S 2 , S 3 , S 2 , S 3 and those satisfying 
Q 2 S [0.7152,0.9] can do it following the sequence S 3 , S 2 , si, S 2 , S 3 . Note that 
from the interval [0.7552, 0.9] both sequences can be taken. 

If not all points in P are “covered” by the fc-length sequences found in the first 
iteration, we restart the procedure from (s, P') where P' C P is the subset of P 
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consisting of the points not covered yet. In our example P' consists of the points 
satisfying q 2 G [0.7, 0.7152]. In the second iteration we find out that the goal orbit 
can be reached from any point in P' by either one of the three 6-transition se- 
quences S3, S2 , S3, S2, S3, S2, S3, S3, S2, S3, S2, Si, S2 , S3 and S3, S2, Si, S2, si, S2, S3, 

and this concludes the computation. The fact that Q 2 does not matter here is 
particular to this example — with other sets of parameters the partition of the 
initial set did involve conditions on q 2 - The reachable states which correspond 
to the discovery of the sequence S3, S2, si, S2, si, S2, S3 in the second iteration are 
depicted in Figure E] and |S1 




Fig. 3. The first iteration of the search tree. The goal orbits were first reached after 
4 transitions along two paths of the tree. 



5 Conclusion 



We have investigated a new methodology for designing hybrid controllers which 
is partially-supported by our reachability analysis tool d/dt. Like |ABD+00j and 
j lTLFOQj this work explores the contribution of the hybrid automaton model to 
the alternative formulation and solution of problems in switching-based control. 
In this paper we have treated an interesting and open problem in robot control 
and provided a partial solution. To improve the performance of the algorithm, we 
plan to investigate other search procedures (backward computation and heuristic 
search) and validate our results via simulation. 
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S3 

Fig. 5. Computation of reachable states for the sequence S3, S2, si, S2, si, S2, S3 contin- 
ued from Figure E] 
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Abstract. In this paper we develop an algorithm for solving the reach- 
ability problem of two-dimensional piece- wise rectangular differential in- 
clusions. Our procedure is not based on the computation of the reach-set 
but rather on the computation of the limit of individual trajectories. A 
key idea is the use of one-dimensional affine Poincare maps for which we 
can easily compute the fixpoints. As a first step, we show that between 
any two points linked by an arbitrary trajectory there always exists a 
trajectory without self-crossings. Thus, solving the reachability problem 
requires considering only those. We prove that, indeed, there are only 
finitely many “qualitative types” of those trajectories. The last step con- 
sists in giving a decision procedure for each of them. These procedures 
are essentially based on the analysis of the limits of extreme trajectories. 
We illustrate our algorithm on a simple model of a swimmer spinning 
around a whirlpool. 



1 Introduction 

One of the main research areas in hybrid systems is reachability analysis. It 
comprises two (closely related) issues, namely, the study of decidability and the 
development of algorithms. Most of the proved decidability results are based 
on the existence of a finite and computable partition of the state space into 
classes of states which are equivalent with respect to reachability. This is the 
case for timed automata 0, and classes of rectangular automata m and hybrid 
automata with linear vector fields uni Except for timed automata, these results 
rely on stringent hypothesis such as the resetting of variables along transitions. 

Although analysis techniques based on the construction of a finite partition 
have been proposed Q, mainly all implemented computational procedures resort 
to (forward or backward) propagation of constraints, typically (unions of convex) 
polyhedra or ellipsoids 111316191111141 . In general, these techniques provide semi- 
decision procedures, that is, if the given final set of states is reachable, they will 
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terminate, otherwise they may fail to. This is a property of the techniques, not 
of the problem, that is, it does not imply that the reachability problem itself is 
undecidable, but only that they do not implement a decision procedure for it. In 
other words, these algorithms may be unsuccessful (i.e., not terminate) for cer- 
tain classes of systems for which the reachability problem is indeed decidable (by 
other means). Nevertheless, they provide tools for computing (approximations 
of) the reach-set for large classes of hybrid systems with linear and non-linear 
vector fields. 

Maybe the major drawback of set-propagation, reach-set approximation pro- 
cedures is that they pay little attention to the geometric properties of the specific 
(class of) systems under analysis. To our knowledge, in the context of hybrid sys- 
tems there are two lines of work in the direction of developing more “geometric” 
approaches. One is based on the existence of (enough) integrals and the ability 
to compute them all fTlll^ . These methods, however, do not necessarily result 
in decision procedures (they are actually not meant to). The other, applica- 
ble to two-dimensional dynamical systems, relies on the topological properties 
of the plane, and explicitly focuses on decidability issues. This approach has 
been proposed in PI- There, it is shown that the reachability problem for two- 
dimensional systems with piece-wise constant derivatives (PCD) is decidable. 
This result has been extended in |5| for planar piece- wise Hamiltonian systems. 
In 0 it has been shown that the reachability problem for PCD is undecidable 
for dimensions higher than two. 

In this paper we develop an algorithm for solving the reachability problem 
of two-dimensional piece-wise rectangular differential inclusions. As in PI , our 
procedure is not based on the computation of the reach-set but rather on the 
computation of the limit of individual trajectories. A key idea is the use of one- 
dimensional affine Poincare maps for which we can easily compute the fixpoints. 
The decidability result of fundamentally relies on the determinism of PCD 
which implies that planar trajectories do not intersect themselves. This property 
is no longer true for differential inclusions. As a first step, we show that between 
any two points linked by an arbitrary trajectory there always exists a trajectory 
without self-crossings. Thus, solving the reachability problem requires consider- 
ing only those. We prove that, indeed, there are only finitely many “qualitative 
types” of those trajectories. The last step consists in giving a decision procedure 
for each of them. These procedures are essentially based on the analysis of the 
limits of extreme trajectories (which do not cut themselves). 

2 Simple Planar Differential Inclusions 

A simple planar differential inclusion system (SPDI) consists of a partition of 
the plane into convex polygonal regions, together with a differential inclusion 
associated with each region. As an example consider the problem of a swimmer 
trying to escape from a whirlpool in a river. 

Example. The dynamics x of the swimmer around the whirlpool is approximated 
by the piece-wise differential inclusion defined as follows. The zone of the river 
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nearby the whirlpool is divided into 8 regions , . . . , ■ To each region Ri 

we associate a pair of vectors (a^jb^) meaning that x belongs to their positive 
hull: ai = bi = (1,5), aa = b 2 = (- 1 , 5 ). ^3 = (“I. in) and bg = (-1,-i), 
&4 = b4 = (-1,-1), as = bg = (0,-1), ag = be = (1,-1), ay = br = (1,0), 
ag = bg = (1, 1). The corresponding SPDI is illustrated in Fig.^ □ 




Fig. 1. The SPDI of the swimmer. 



More formally, a SPDI is a pair TL = {V, 4>), where P is a finite partition of 
the plane into convex polyhedral sets, and for each P G P, 4>{P), also denoted 
by is the set of all linear combinations x = a ap + f3 bp, with a,f3 > 0, 
and a + /3 > 0, of two vectors ap and bp, such that ap • bp < 0, where • is the 
scalar product and ap = ( 02 ,- 01 ) is the clockwise rotation of ap by the angle 
^ (notice that ap • ap = 0). 

Let E{P) be the set of edges of P, that is, the set of open segments forming 
the boundary of P, and V{P) be the set of vertices in the boundary of P. We 
say that e is an entry of P if for all x S e and for all c S ^(P), x + ce S P 
for some e > 0. We say that e is an exit of P if the same condition holds for 
some e < 0. We denote by in(P) C E{P) the set of all entries of P and by 
out(P) C E{P) the set of all exits of P. In general, E{P) ^ in{P) U out(P). 
We say that P is a good region iff all the edges in E{P) are entries or exits, 
that is, E{P) = in{P) U out(P). Notice that, if P is a good region, then for all 
e G E{P), the director vector of e does not belongs to 4>{P) (Fig. 0. Hereinafter, 
we assume that all regions are good regions. Let x G V(P) be a common vertex 
of two edges e and e'. x is an entry point to P if both e and e' are entry edges; 
it is an exit point if both e and e' are exit edges. In fact, vertices can be seen as 
a particular kind of edges, with exactly one point. In what follows the term edge 
will be understood as belonging to the set EV{P) = E{P) U V{P). If needed, 
the difference between edge and vertex will be explicitly specified. 

A trajectory in some interval [0,P] C K, with initial condition x = xq, is a 
continuous and almost-everywhere (everywhere except on finitely many points) 
derivable function ^(-) such that ^(0) = xg and for all t G (0,P), if ^{t) G 
P\ EV{P), then ^{t) is defined and ^{t) G 4>{P)- 

The point-to-point reachability problem for R, is the following: Given x, x' G 
K^, is there a trajectory ^ and t > 0 such that ^(0) = x and ^{t) = x'?. If the 
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Fig. 2. a) A good region, b) A bad region (es ^ in{P) U out{P)). 



answer is yes, we say that x' is reachable from x. The edge-to-edge reachability 
problem is the following: Given two edges e and e' oi PL, is there x S e and x' G e! 
such that x' is reachable from x? The region-to-region reachability problem is 
defined similarly. 



3 Properties of Trajectories 

W.l.o.g. we will consider in what follows that ^(0) G e for some edge e. The 
trace of a trajectory ^ is the sequence r(^) = XgXi ... of the intersection points 
of ^ with the set of edges, that is, x^ G ^ fl [}EV{P) for all P GV. The edge 
signature of ^ is the sequence (t(^) = eoei . . . of traversed edges, that is, Xi € Ci. 
The region signature of ^ is the sequence p(^) = PqPi ... of traversed regions, 
that is, Ci G in{Pi). 

Let ^ be a trajectory whose trace is r(^) = Xq...x/c. Let 0 = tp < G < 

. . . < tfc be such that ^(ti) = Xj. Since ^ is continuous and derivable in the 
interval (ti,ti^i), there exists a unique trajectory with ^'(ti) = ^(ti) for all 
i G [0, k— 1], such that the derivative is constant in the interval {ti,ti+i). That 
is. 

Proposition 1. For every trajectory ^ there exists a trajectory with the same 
initial and final points, and edge and region signatures, such that for each Pi 
in the region signature, there exists Ci G 4>{Pi), such that f'(t) = for all 
t S {ti, ti+i). 

Hence, in order to solve the reachability problem it is enough to consider trajec- 
tories having piecewise constant slopes. Notice that, however, such slopes need 
not be the same for each occurrence of the same region in the region signature. 
Hereinafter, we use the word “trajectory” to mean trajectories whose derivatives 
are piecewise constant. 

Consider a region P and let c G (j>{P)- The mapping 17 : K., defined 

as f7(x) = X • c, assigns to every x G a value proportional to the length of 
the projection of the vector x on the right rotation of c (see 0). Indeed, the 
ordering is given by the direction of c and one can easily see that the relation 

defined as xi A X 2 if l7(xi) < f?(x 2 ), is a dense linear order on in{P) and 
out{P) (Fig. EJ. We use ^ to denote the strict variant of ^ and say that ei ^ 62 
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iff Cl ^ 62 and Xi ^ X2 for every Xi G ei,X2 G 62- For example, in Fig. 0 we 
have 6i ^ 62 ^ 63. Notice that the order does not depend on the choice of c. 




Fig. 3. Ordering: xi ^ X 2 . 



We say that a trajectory ^ crosses itself if there exist such that ^{t) = 

If a trajectory does not cross itself, the sequence of consecutive intersection 
points with in(P) or out(P) is monotone with respect to That is, for every 
three points xi, X2 and X3 (visited in this order), if xi ^ X2 ^ X3 the trajectory 
is a “counterclockwise expanding spiral” (Fig. 0 (a)) or a “clockwise contracting 
spiral” (Fig. 0 b)) and if X3 ^ X2 ^ xi, the trajectory is a “counterclockwise 
contracting spiral” (Fig. 0 c)) or a “clockwise expanding spiral” (Fig. 0 d)). On 
the other hand, if the sequence of intersections points with in{P) or out{P) is 
monotone (both increasing or both decreasing), the trajectory does not cross 
itself. 

Lemma 1 . For every trajectory if ^ does not cross itself, then for every edge 

e, the sequence r(^) C\e is monotone (with respect to ^). 





Fig. 4. Spirals. 

Now suppose that the trajectory f with trace T{f) = xq . . . xy crosses itself once 
inside the region P. Let 61, 62 G in{P) be the input edges and e[,e'2 G out{P) be 
the output ones. Let x = Xj G 61 and y = Xj G 62, with i < j, be the points in 
r(^) the first and the second times f enters P, and let x' = x^+i G 63 and y' = 
Xj+i G e[ be the corresponding output points. Let Cx,Cy G (j>{P) = be the 
derivatives of f in the time intervals (ti, ti+i) and (tj,tj+i), respectively. Indeed, 
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Cx and Cy are the director vectors of the segments xx' and yy', respectively 
(Fig. Eta)). 

Consider now the segment xy'. Notice that the director vector of this 
segment can be obtained as a positive combination of the vectors and Cy. 
Thus, S Hence, there exists a trajectory that does not cross itself in 

P having a trace r(^') = xq . . . xy' . . . x/ (Fig. Elb)). Notice that the result also 
works for the degenerate case when the trajectory crosses itself at an edge (or 
vertex) . 




Fig. 5. Obtaining a non-crossing trajectory 



If the trajectory ^ crosses itself more than once in region P, then the number 
of times the trajectory obtained by cutting away the loop (Fig. Etc)), crosses 
itself in P is strictly smaller than the number of times ^ does it (see Fig. 0. 
After replacing xx' and yy' by xy', the intersection q of xx' and yy' disappears. 
If the new segment of line xy' crosses another segment zz' (say at a point t), 
then zz' necessarily crosses either xx' (at r) or yy' (at s) -or both-, before the 
transformation. The above is due to the fact that if zz' crosses one side of the 
triangle xy'q then it must also cross one of the other sides of the triangle, say 
at r. Thus, no new crossing can appear and the number of crossings in the new 
configuration is always less than in the old one. 





Fig. 6. The number of crossings decreases: (a) Before (3 crossings); (b) After (1 cross- 
ing). 
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Lemma 2. For every trajectory ^ that crosses itself at least once, there exists 
a trajectory with the same initial and final points of ^ having a number of 
self- crossings strictly smaller. 

The above result follows from a straightforward inductive reasoning, as well as 
the following one. 

Proposition 2. If there exists an arbitrary trajectory from points Xq G Cq to 
Xf € 6f then there always exists a non-crossing trajectory between them. 

Hence, in order to solve the reachability problem we only need to consider non- 
crossing trajectories with piecewise constant derivatives. In what follows, we only 
deal with trajectories of this kind. 

4 Properties of Edge Signatures 

Let ^ be a trajectory with trace r(^) = xq . . . Xp, edge signature (t(^) = cq . . . Cp, 
and region signature p(f) = Pq . . . Pp. An edge e is said to be abandoned by 
^ after position i, if = e and for some j, k, i < j < k, Pj . . . Pk forms a 
region cycle and e ^ {e^+i, . . . , e^}. Since trajectories are finite we should add 
the trivial case when e yf Cj for all j > i. 

Lemma 3 (Claim 2 in |3)- Aor every trajectory ^ and edge e, if e is abandoned 
by f after position i, e will not appear in tT(^) at any position j > i. 

Given a sequence s, we use notations first{s) and last{s) for the first and last 
elements of the sequence respectively, e denotes the empty sequence An edge 
signature cr(^) can be canonically expressed as a sequence of edges and cycles of 
the form = risJV 2 S 2 ^ . . . where 

1. For alH S [1, n -b 1], is a sequence of pairwise different edges; 

2. For all i S [l,u], Si is a simple cycle (i.e., without repetition of edges) 
repeated ki times; 

3. For all i € [l,n — 1], first(ri+i) yf first{si) if yf e, otherwise 

first{si+i) yf first(si); 

4. For all i G [l,n], if yf e then last{ri) = last{si); 

5. r„_|_i yf e. Moreover, Vn+i = first(sn) if cr(^) ends in a loop and 
first{rn-i-i) yf first{sn) otherwise. 

This canonical representation can be obtained as follows. Let cr(^) = e\ . . . Cp-iCp 
be an edge signature. Starting from Cp-i and traversing backwards, take the 
first edge that occurs the second time. If there is no such edge, then trivially 
the signature can be expressed in a canonical way and we are done. Otherwise, 
suppose that the edge Cj occurs again at position i (i.e. = Cj with i < j), 

thus (Tc(^) = wsr, where w,s and r are obtained as follows, depending on the 
repeated edge: 

w = eo. . .6i 

S — . . . €j 

T — . . . 6p_i 
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Clearly r is not a cycle and s is a simple cycle with no repeated edges. We 
continue the analysis with w. Let km = max {I \ sMs a suffix of ru}. Thus, 
(Tc(C) = w's^r with w' = . . .Ch (a prefix of w) and k = km + 1- We repeat 

recursively the procedure above with w' . Adding the edge Cp to the last r (at the 
end) we obtain ctc(^) = . . .r„s^"r„+i that is a canonical representation of 

signature a. 

Let us define the type of a signature a as type{a{f)) = ri, si, . . . , r„, s„, r„+i. 
Notice that the “preprocessing” (taking away the last edge Cp) is done in order 
to differentiate edges signatures that end with a cycle from those that do not. 
There exists many other (maybe easier) ways of decomposing a signature cr in a 
canonical form (in particular, traversing forward instead of backwards), but the 
one chosen here permits a clearer and simpler presentation of the reachability 
algorithm. In fact in this canonical form, the last visited edge in a cycle ei . . . 
is always the last one (cfc). 

Example. Let us consider the following examples. Suppose that a = abcdbcefg 
efgefgefhi. Then, after applying once the above procedure we obtain that CTc = 
ic(s 2 )^?’i, with w = abcdbcef; S 2 = gef; ri = h. Applying the procedure once 
more to w we obtain w = w'{ss)^r 2 with w' = = abc] S 3 = dbc; T 2 = ef. 

Putting all together and adding the last edge (t) gives CTc = abc{dbc)^ef{gef)^hi 
with type type{<j) = abc, dbc, ef, gef, hi. Suppose now, that the signature ends 
with a cycle: a = abcdbcef gef gef gef gef . In this case we apply the preprocessing 
obtaining CTc = w{s 2 )'^ri with w = abcdbce] S 2 = fge; r\ = £. Applying the 
procedure to w we finally obtain w = w'{s 3 )^r 2 with w' = = abc; S 3 = 

dbc; V 2 = e and that gives CTc = abc{dbcye{fge)‘^f (adding / to the end). □ 

Lemma 4. The type of a signature a, type{a), has the following properties: 

1. For every 1 < i yf j < n + 1, and rj are disjoint; 

2. For every ^ < i ^ j < n, Si and sj are different; 

3. If V is a vertex appearing in type{a), then it can only occur exactly once in 
Vi for some l<i<n+l in a. Moreover, v ^ last{ri) unless i = n+1. 



Proposition 3. The set of types of edge signatures is finite. 

Thus, to solve the reachability problem we can proceed by examining one by one 
the types of signatures. 

5 AfRne Operators 

Before getting into the problem of analyzing types of edge signatures, we need 
to introduce some useful notions. 

An affine function / : R — >■ R is defined by a formula f{x) = ax + b with 
a > 0. An affine multivalued operator F : R — >■ 2® is determined by two affine 
functions fi{x) and fu{x) and maps x to the interval {fi{x), /„(x)), where (a, b) 
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means (a,b), [a,b], {a,b] or [a, 6) : F{x) = We use the nota- 

tion F = (fi,fu)- Such an operator can be naturally extended to subsets of 
R: F{S) = Ua:GS^(^)- particular, if S' = {l,u): F{{l,u)) = (/;(0. /«(«))• A 
truncated affine multi-valued operator G : R — >■ 2* is determined by an affine 
multi-valued operator F and an interval {L, U) as follows: G{x) = F(a;) fl (L, U). 
Such operators can be also extended to sets. We use notations G = F (1 {L,U) 
and F = G. 

Lemma 5 (composition of afRne operations). Affine functions, affine 
multi-valued operators, and truncated affine multi-valued operators are closed 
under composition. 

Example. Let Gi(a:) = (2cc-|-3, 3x-|-5] and G 2 {x) = [bx-\-2, 7a; -I- 6 ] be two (non- 
truncated) affine multi-valued functions, Gi = Gin(l, 6 ], and G 2 = G 2 n[ 6 , 10 ) 
their truncated versions. The truncated affine multi-valued operator G 2 o G\{x) 
is obtained as follows: 

G 2 o Gi(x) = G 2 o Gi{x) n G 2 (( 1 , 6 ]) n [ 6 , 10) 

= (5(2a: -b 3) -b 2, 7(3a: -b 5) -b 6] n (5 • 1 -b 2, 7 • 6 -b 6] n [6, 10) 

= (10a; -b 17, 21a; -b 41] n (7, 48] n [6, 10) 

= (10a; -b 17, 21a; -b 41] n (7, 10). 

Notice also that for any interval {l,u) its image is G 2 o Gi{{l,u)) = (10/ -b 
17, 21u-b41) n(7, 10). □ 

Let / be an affine function, xq be any initial point and a;„ = /"(a;). Clearly, 
the sequence a;„ is monotonous, and it converges to a limit x* (finite or infinite) . 
Indeed, x* can be effectively computed knowing a, b and xq, as follows. If a < 1, 
a;* is the unique fixpoint of /, that is, ax* -hb = x*, which yields x* = 6/(1 — a). 

If a = 1, a;* = —00 if 6 < 0, a;* = 00 if 6 > 0, and x* = xq, if 6 = 0. If a > 1, 

let a;* = 6/(1 — a), then x* = —00 if a;o < a;*, a;* = 00 if xq > a*, a;* = a;o = a;*, 
otherwise. This result can be extended to multi-valued affine functions. 

Lemma 6. Let (lo,Uo) be any initial interval and {ln,Un) = F”‘{{Iq,uq)). Then 

1. The sequences In and are monotonous; 

2. They converge to limits I* and u* (finite or infinite ), which can be effectively 
computed. 

Proposition 4. Let F be truncated affine and L C {L,U). Then F"‘{I) = 
F-{L) fML,U). 

6 Computing the Successor Function 

To solve the reachability problem for SPDI, the next step is to provide a pro- 
cedure for computing the successors of a point (and an interval), which requires 
having an effective representation of (rational) points and intervals on edges. 
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Let us first introduce a one-dimensional coordinate system on each edge. For 
each edge e we chose (1) a point on it (the origin) with radius-vector v, and 
(2) a director vector e going in the positive direction in the sense of the order 
Now to characterize e we need the coordinates of its extreme points, namely, 
e^e“ G QU{— 00 , 00 } such that e = {x = v-|-xe | < a; < e“|. That is, an edge 

e £ E can be represented by a triplet (v, e, (e*, e“)). If the edge is a vertex, the 
representation is simply (v, [0, 0]). Now, every point x = v-|-xe G e is represented 
by the pair (e,x) (Fig[3(a)), and every interval (xi,X 2 ) C e is represented as 
(e, (xi, X 2 )), where Xi = (e,Xi) and X 2 = (e,a; 2 ) (Figl^b)). Now, having fixed 





Fig. 7. (a) Representation of edges; (b) Representation of an interval; (c) One-step 
successor. 



a one-dimensional coordinate system to represent points, the question now is to 
take advantage of it to compute the successor of a point or an interval. 

Let e = (e^e“) G in{P) and e' = (e'\e'“) G out{P). For x = (e,x) and 
c G (/i(F’), we denote by SucCg_g,(a;) the unique x' = (e', x') such that x' = x-1- cf 
for some t > 0. The point {e',x') is the successor of (e,x) in the direction c (see 
FigQc)). Expanding, v'-|-x'e' = v-bxe-l-te. Multiplying both expressions by c 
we obtain that (v'-l-x'e')c = v-c-|-a;ec, i.e. x'{e' -c) = x(e-c)-|-(v — v')-c. Thus, 
x' = SucCg g/(a;) = -c and x' G (e'^ e'“). Indeed, putting o(c) = 

and /3(c) = ■ c we have the following result. 

Lemma 7. The function SucCgg/(a;) = a{c)x + /3(c) n is truncated 

affine. 

SucCg^g'(a;) will denote the non-truncated function a(c)x + /3(c). The notion 
of successor can be extended on all possible directions c G 4>(P) and it can be 
applied to any subset S C (e*,e“) and in particular to intervals (l,u): 

Lemma 8. Let 4>{P) = x = {e,x) and {l,u) C (e^e“). Then: 

1 . Succg,g/(a;) = Ucg0(p) SucCg_g,(a;) = (SucCg_g,(x), SucCg_g,(a;)) n (e'\ e'“); 

2. SucCg, £/((/,«)) = (SucCg_g,(/), SucCg_g/(u)) n (e'\ e'“). 

The successor operator will be used as a building block in the reachability al- 
gorithm. It can be naturally extended on edge signatures: for w = ei, 62 , . . . , e„ 
let 

Succ„,(/) = Succg„_i,g„ o . . . o Succg2,g3 o SucCg,,g2(J) 
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that by Lemma 0 is truncated affine. Notice that since we use edge signatures 
the semi-group property takes the following form. 

Lemma 9. For any edge signatures w and v and an edge e, SucCeu, o Succ„e = 
Succ^g^ . 



Example. Let us come back to the example of the swimmer trying to escape 
from a whirlpool in a river (see Fig. Suppose that the swimmer is following 
a trajectory with edge signature (ci . . .eg)*. It is not difficult to find a repre- 
sentation of the edges such that for each edge e,, (e(,e“) = (0, 1). Besides, the 
truncated affine successor functions are: 



[X x~\ 

Succg,e3(a:) 

SucCgjei+i (®) = [®, x] n (0, 1), for all i e [3, 7] SucCggei (®) 



3 2 

*-f0’*+15 
1 1 

x+ -,x+ - 
5 5 



n(o,i) 
n(o,i) 



The successor function for the loop s = ei . . . eg is obtained by composition of 
the above functions as follows. Let us first compute 



SuCCg3g2g3 (^, w) SuCCg2g3 O SuCCg3g2 



(l,u) 



- [5 “ I + 1 + A) (O’ 1) 



10 ’ 

10 ’ 



15J 

A' 

15 



n(o,i) 



Since Succ 



EiEi+l 



for i S [3, 7] are the identity functions, we have that 



SucCp 



gg(/,'u) SuCCggg^ O SuCCe3e2g3 



(l,u) 



= \- — ^-i-i - i-i-ilnlo 11 



10 

1 u 
10 ’ 2 






By Lemma El we have that I* 



- and u* = jfx 



7 Reachability Analysis 



The algorithm for solving the reachability problem between two points xq = 
(eo,a:o) and xy = (e/,x/) is depicted in Fig. 0 The proofs of soundness and 
termination are given in the extended version (|S|). It works as follows. 



Reach. From the section above we know that there exists a finite number of type 
signatures of the form ri, si, . . . , r„, s„, Xn+i. Moreover, the type signatures are 
restricted to those with Cq = first{ri) and ey = last{rn+i) . Given such a set of 
type signatures type(eo, ey), the algorithm Reach{-) is guaranteed to terminate, 
answering YES if xy is reachable from Xq or NO otherwise. Reachability from 
xo to Xy with fixed signature w is tested by the function Reachtype{xo,Xf,w). 
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Reachtype- Let the type w have the form w = ri, si, . . . , r„, s„, rn+i. Put 
fi = first{si) and exi = first{ri+i) if is non-empty and /i+i otherwise 
(i.e. eXi is the edge to which the trajectory exits from the loop Si). Let us say 
that a type signature w has a loopend property if r„_|_i = first(sn), i.e. signa- 
tures of type w terminate by several repetitions of the last loop. This algorithm 
uses two functions: Test{S,s,x) that answers whether x is reachable from a set 
S (represented as a finite union of intervals) in the loop s (formally, whether 
X G SucCg+ first{s){I))', and the function Exit{S, s, e) that for an initial set S, a 
loop s, and an edge e (not in this loop) finds all the points on e reachable by 
making s several times and then exiting to e (formally, it computes SucCs+e(L), 
which is always a finite union of intervals ) . Since we know how to calculate the 
successor of a given interval in one and in several steps (SucCee'(-) and SucCr(-))i 
in order to implement Test(-) and Exit{-) it remains to show how to analyze 
the (simple) cycles Si and eventually their continuation. Both algorithms Test(-) 
and Exit{-) start by qualitative analysis of the cycle implemented in the function 
Analyze{I,s). This analysis proceeds as follows. 

Analyze. The function Analyze{I, s) returns the kind of qualitative behavior of 
the interval I = (l,u) under the loop s. Let s be a simple cycle, / = first(s) its 
first edge, and I = {l,u) C f an initial interval and SucCsj(a;) = SucCs,/(x) H 
{L,U). The first thing to do is to determine the qualitative behavior of the 
leftmost and rightmost trajectories of the interval endpoints in the cycle. This 
can be done without iterating SucCg/. Indeed, by LemmaCJQ we can compute the 

n 

limits = limji_).oo Succ^ j((^, it)) (notice that those are limits only for the 

non-truncated operator Succ), not taking into account that the edges are possible 
bounded (we use Lemma 0 and compare these limit points corresponding to 
unrestricted dynamics with L and U. There are five possibilities: 

1. STAY The cycle is not abandoned by any of the two trajectories: L <1* < 
u* < U. 



function fteac/i(xo, x/) 

R — false 

for each w E type{eo,ef) 

R = R\/ Reachtype{xo,Xf,w) 
i — R 


function Reachtype{xo,Xf,w) : 

S = SuCCrifi (xo) 
for i = 1 to n — 1 

S = SucCn^ifi^j^{Exit{S, Si, exi)) 

if loopend(u') 

then i — Teat{S,Sn,Xf) 

else i Xf E SuCCr„^,^{Exit{S,Sn,eXn))'l 


fnnction Exit{S, s, ex) 

F = 0 

for each I E S such that SucCg /(/) / 0 
E — £/ U ExitAnalyze{^^eCs /(f), Cx) 
< — E 


function Test{S, s,x) 

R = false 

for each I E S such that SucCs f{I) 

R — R \/ T est Analyze (SuCCs j (f ) , S, x) 
< — R 



Fig. 8. Main algorithm. 
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2. DIE The right trajectory exits the cycle through the left (consequently the 

left one also exits) or the left trajectory exits the cycle through the right 
(consequently the right one also exits). In symbols, u* < L\/ I* > U . 

3. EXIT-BOTH Both trajectories exit the cycle (the left one through the left 

and the right one through the right ) : I* < L Au* > U. 

4. EXIT-LEFT The leftmost trajectory exits the cycle but not the other: I* < 

L<u* <U. 

5. EXIT-RIGHT The rightmost trajectory exits the cycle but not the other: 

L<1* <U <u*. 

Exit. The function Exit{S, s,ex) should return SucCg+ea;(5'). Both the argu- 
ment S and the result are finite collections of intervals. The exploration is made 
for each initial interval separately. Notice that the call SucCg j{I) ensures that 
I C {L,U). All the work for each initial interval I is done by the function 
ExitAnaiyze{lTS,ex) which launches the Analyze{-) procedure described above 
and last, according to the result of this analysis launches one of five special- 
ized procedures ExUstay , Exitj^EFT, ExUhight, ExUboth, ExUjjie which 
calculates the exit set (Fig. dj. 



function Found{I, x) 



function Search{I, x) 
while Found{I, x) = NOTYET 
I = SuCCs /(/) 
i — Founded, x) 



cases 
X € I : 

/ = 0 : 
a: < / A f t : 

X > I Au i : 
else : 




YES 

NO 

NO 

NO 

NOTYET 



Fig. 9. Search and Found. 



Test. The upper-level structure is the same as for EXIT: each initial interval 
is treated separately by Test Analyze, which makes one turn of the loop, calls 
Analyze and delegates all the remaining to one of the five specialized functions 
TestsTAY, TestLEFT, TestRiGHT, TestsoTH, TestniE (Fig. dj). The five spe- 
cialized Test functions use the following two procedures (Fig. EJ. The function 
Eound{I, x) determines if the current interval I contains x (YES), does not con- 
tain X and moves in the opposite direction (NO), or none of both these cases 
(NOTYET). Found{I,x) uses the fact that the sequences In and are increas- 
ing or decreasing (which can be easily determined at the stage of the preliminary 
analysis of the loop): I f means that the sequence l,li,l 2 , - ■ ■ of successive suc- 
cessors of I is increasing whereas I j, means that the sequence is decreasing, and 
similarly for u f and u j,. The function Searched, x) iterates the loop s until the 
previous function Found gives a definite answer YES or NO (Fig. ED . It is used 
only when its convergence is guaranteed. 
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Exit 


Test 


STAY 


function Exit stay {I ^s, ex) 
^ 0 


function TestsTAvil, s,x) 

cases 

1* <x<u* :i — YES 
X < T A 1 X ’. < — NO 
X > u* A uX '■< — NO 
else : < — Search{I, x) 


DIE 


fnnction ExitoiE{I,s,ex) 
f = first{s) 

S = 0 

repeat 
I = SuCCs/(l) 
S^S-USucc. , 41 ) 
until 1 = 0 
< — S 


fnnction TestoiE{I,s,x) 
i — Search{l, x) 


BOTH 


function ExitsoTHi.! , s,ex) 
i SuCCs ea:((L, U)) 


function TestBOTH{I,s,x) 
^xe {L,U)7 


LEFT 


function ExitLEFT{I , s,ex) 
i — Succs ,x{{L,u)) 


function TestLEFril , s, x) 

cases 

x£{L,u*)-. < — YES 

x<{L,u*): i — NO 

(L,u*) < X AuX : — NO 
else : < — Search(I, x) 


RIGHT 


Similar to the previous case. 


Similar to the previous case. 



Fig. 10. Exit and Test. 




Fig. 11. Example: x/ = (ei, |) is not reachable from xo = (ei, |) {u* < |). 



Example. Consider again the swimmer. Let xq = (ei, be her initial position. 
We want to decide whether she is able to escape from the whirlpool and reach 
the final position x/ = (ei, |). Recall that I* = | and u* = = |. 

Thus, by the Analyze function we know that the cycle behaves as an Exit-LEFT 
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and applying the function TestLEFT we obtain that x/ = (ei, |) is not reachable 
from xo = (ei, because we have that u f and u* < Xf (| < |) (Fig. El- D 
From all the results above we have the following theorem. 

Theorem 1 (Point— to— Point Reachability). The point-to-point, edge-to- 
edge and region-to-region reachability problems for SPDI systems are decidable. 

□ 



8 Concluding Remarks 

We have presented an algorithm for solving the reachability problem for simple 
planar differential inclusions. The novelty of the approach for the domain of Hy- 
brid System is the combination of two techniques, namely, the representation of 
the two-dimensional continuous dynamics as a one-dimensional discrete system 
(due to Poincare), and the characterization of the set of qualitative behaviors of 
the latter as a finite set of types of signatures. 

One possible direction of future work is to try to apply the same method for 
solving the parameter synthesis problem for SPDFs, that is, for any two points, 
xq and xy, assign a constant slope cp € 4>{P) to every region P such that xy is 
reachable from xq, or conclude that such an assignment does not exist. Clearly, 
the decidability of the reachability problem does not imply the decidability of 
the parameter synthesis one. 

Another question that naturally arises is decidability of the reachability prob- 
lem for hybrid automata whose locations are equipped with SPDFs. We can cer- 
tainly find (stringent) conditions, such as planarity of the automaton, “memory- 
less” resets, etc., under which decidability follows almost straightforwardly from 
the decidability of SPDFs. On the other hand, it is not difficult to see that 
this problem, without such conditions, is equivalent to deciding whether given 
a piece-wise linear map / on the unit interval and a point x in this interval, 
the sequence of iterates x, f{x), f{f{x)), and so on, reaches some point y. This 
last question is still open ng. And last but not the least, another interesting 
issue is the complexity analysis of the algorithm. It should be based on counting 
all “feasible” types of signatures. Our finiteness argument of lemma 0 gives a 
doubly exponential estimation, which can certainly be improved. 
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Abstract. The behavior of the run of an impulse differential inclnsion, 
and, in particular, of a hybrid control system, is “summarized” by the 
“ initialization map” associating with each initial condition the set of 
new initialized conditions and more generally, by its “substratum” , that 
is a set-valued map associating with a cadence and a state the next 
reinitialized state. These maps are characterized in several ways, and 
in particular, as “set-valued” solutions of a system of Hamilton-Jacobi 
partial differential inclusions, that play the same role than usual 
Hamilton-Jacobi-Bellman equations in optimal control. 

Keywords: hybrid control, impulse control, differential inclusion, vi- 
ability, run, execution, periodic, cadenced run, equilibrium, Kakutani 
Theorem, contingent cone, Marchaud map. 



1 Introduction 

Impulse differential inclusions, and in particular, hybrid control systems, are de- 
fined by a differential inclusion (or a control system) and a reset map. A run of an 
impulse differential inclusion is defined by a sequence of cadences, of reinitialized 
states and of motives describing the evolution along a given cadence between two 
distinct consecutive impulse times, the value of a motive at the end of a cadence 
being reset as the next reinitialized state of the next cadence. 

A first advantage of introducing impulse differential inclusions is to sum- 
marize the usually protracted description of an hybrid systeirQ by only two 
set-valued maps F — the right-hand side of the differential inclusion governing 
the continuous evolution of a hybrid system — and i?, describing the reset map 
reinitializing the system when required and a constrained set K inside which 
the evolution of the “run” or “execution” must remain. Hence, for instance, the 
existence of a run of an hybrid system for every initial set becomes a viability 
problem of an adequate auxiliary subset under an impulse differential inclusion, 
that can be characterized elegantly end efficaciously. 

^ See for instance among many papers and books |1 31 Branicky, Borkar & Mitter], |1 21 
Bensoussan & Menaldi], |'l 711 iSl Matveev & Savkin] and |2()l Shaft & Schumacher]. 
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The behavior of the run is “summarized” by the “ initialization map” U := 
U(F,R) associating with each initial condition xq € K the set of new initialized 
conditions Xi £ B(x(~ti)) when x(-) ranges over the set of solutions to the 
differential inclusion x' € F{x) viable in K until they reach at time 

> 0 at x{~ti) G R~^{K). 

Indeed, the sequence of successive initial conditions of a viable run x(-) 
of the impulse differential inclusion {F, R) — constituting the “discrete compo- 
nent of the run” — is governed by the discrete system G U(^F,R){xn-i) H K 
starting at Xq. The knowledge of the sequence of initialized states x„ allows us 
to reconstitute the “continuous component” of the run by solving the differen- 
tial inclusion x' G F{x) starting at each reinitialized state and satisfying the 
end-point condition Xn+i G R{x{~tn+i)), which exists thanks to the definition 
of the map . 

Assume for a while that the impulse differential inclusion is actually an im- 
pulse differential equation (/, r) where the maps / and r are single- valued and 
that the initialization map is single- valued and differentiable. Then we shall 
prove that the initialization map is a solution to the system of first-order partial 
differential inclusions 



V i = 1, . . . n, 



E 



t=i 



duj{x) 

dxj 






= 0 



or, in a more compact form, 0 = — 

ox 

V X G AT n r~^ 



f{x), satisfying the “condition 
(AT), r(x) = u{x) 



5? 



Actually, we shall extend this result to general impulse differential inclusions by 
characterizing the initialization map C^(F,fl) a generalized (set- valued) solu- 
tion — a Frankowska solution — to the system of first-order partial differential 
inclusions 

0 gfW 

satisfying the “condition” 

y X £ KnR-\K), R{x) C U{x) 



These are indeed really Dirichlet boundary condition whenever the reset map 
R is defined only on the boundary dK of a closed subset K and maps dK into 
the interior of K. In this case, resetting initial conditions happens only when 
the continuous evolution of the state governed by the differential inclusion or 
the control system is about to leave the domain K. Hence the reset map assigns 
new initialized states in the interior of K. 

We shall introduce more generally another set- valued map summarizing the 
behavior of an impulse differential inclusion, called the substratum, that is the 
topic of this paper. 
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Outline: We begin by giving our definition of impulse differential inclusions. 
We then recall the characterization of viable subsets under an impulse differential 
inclusion an derive from it a necessary and sufficient condition for the existence 
of solutions to hybrid differential inclusions. Then, we devote the next section to 
the graphical properties of the initialization map U and we derive its properties 
from the general properties of viable-capture basins of a target by a differential 
inclusion. 

In the last section, we translate the Frankowska characterization of viable- 
capture basins in terms of kinds of systems of first-order Hamilton- Jacobi partial 
differential equations characterizing the substratum the solutions of which are 
the initialization maps and the substratum. 

2 Impulse Differential Inclusions 

“Hybrid control systems”, as they are called by engineers, or “multiple-phase 
dynamical economies”, as they are called by economists (see for instance (El 
Day]), or “Integrate and Fire” models in neurobiology (see for instance (l 41 
Brette]) — may be regarded as impulse differential inclusions. 

Here, X := R" and Y := R™ denote finite dimensional vector spaces. Let 
f : X X Y I— >-Jfbea single-valued map describing the dynamics of a control 
system and P : X Y the set-valued map describing the state-dependent 
constraints on the controls. 

First, any solution to a control system with state-dependent constraints on 
the controls 

(i) x'(t) = f{x(t),u{t)) 

\ii) u{t) € P{x{t)) 

can be regarded as a solution to the differential inclusion x'{t) G F{x{t)) where 
the right hand side is defined by F{x) := f{x,P{x)) := {f{x,u)}ueP(x)- 

Therefore, from now on, as long as we do not need to implicate explicitly the 
controls in our study, we shall replace control problems by differential inclusions. 

We shall say that K is locally viable under F if from every x G K starts a 
solution a;(-) to the differential inclusion x' G F(x) viable in K on the nonempty 
interval in the sense 



Vf G [0,T,[, x{t) G K 

and that K is viable if we can take = -|-oo. It is loeally backward invariant 
under F if for every to G]0, +oo[, x G K, for all solutions a:(-) to the differential 
inclusion x' G F{x) arriving at x at time to, there exists s G [0, to[ such that x(-) 
is viable in K on the interval [s,to], and backward invariant if we can take s = 0. 
We denote by 



Graph(F) := {(x,y) G X G Y j y G F(x)} 



the graph of a set-valued map F : X Y and Dom(F) := {x G X\F{x) yf 0} 
its domain. 
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Let us set x{~ t) := x{r) when x{-) is defined on some interval [t—r], t[ 

where rj > 0, and, for consistency purposes, a;(s) = x{~ t) if s = t. An impulse 
differential inclusion (and in particular, an impulse control system) is described by 
a pair (F,R), where the set- valued map F : X X mapping the state space 
X := R" to itself governs the continuous evolution of the system in K and where 
R, the reset map, governs the discrete switches to new “initial conditions” when 
the continuous evolution is doomed to leave K. 

Such a hybrid evolution, mixing continuous evolution “punctuated” by dis- 
continuous impulses at impulse times is called in the “hybrid system” literature 
a "run" or an "execution". 

Definition 21 Let us consider a finite dimensional vector space X, a closed 
subset K G X , a set-valued map F \ X X and a set-valued map R : X X , 
regarded as a reset map. We regard the pair (F, R) as the dynamics of an impulse 
differential inclusion. 

A run of the impulse differential inclusion is a map x(-) from [0,T] to X 
if T < -boo or from [0,-|-oo[ to X if T = -|-oo which is associated with a non 
decreasing sequence T{x{-)) := {tn}n>o o/ impulse or switching times to := 0 < 
ti < ■ ■ ■ < tn < ■ ■ ■ < T (depending on the run x{-) ) such that 

1. xft-YiJ^fj G R(^x(^tji^^ if tn-t-l — tji, 

2. or else, on the interval [f„, tn-i-i[, xf) is a solution to the differential inclusion 
x' G F{x) starting at x(fn) at time tn until time tn+i at which we take 
x{tn+i) G R{x{ ^n+l)) ■ 

We denote by Tn := t„ — tn-i the nth cadence of the run and by x„(-) := x{--\-tn) 
the nth motive of the run, a solution to the differential inclusion x' G F{x) 
starting at x(fn) on the interval [0,t„] and satisfying the end-point condition 
Xn(xn) G Ji~^(Xn-i-i)- The scquencc of states x{tn) is called the sequence o/ ini- 
tialized states. 

We say that a run x{-) is viable in K if for any t > 0, x{t) G K. 

At this stage, a run a:(-) can just be a (discrete) sequence of states Xn+i G 
R{xn) at a fixed time, or just a (continuous) solution a:(-) to the differential 
inclusion x' G F(x), or an hybrid of these two modes, the discrete and the 
continuous. 

Hybrid systems can be regarded as instances of viable impulse differential 
inclusions: we refer to 0 Aubin] or m Aubin, Lygeros, Quincampoix, Sastry 
& Seube] for more details on that topic. 

3 The Substratum and the Initialization and Impnlse 
Maps 

We denote by S{x) C C(0,oo; A) the set of absolutely continuous functions t H> 
x{t) G X satisfying 



for almost all f > 0, x'{t) G F{x{t)) 
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starting at time 0 at x: x(0) = x and by : K C(0,oo,i4T) the viable 
solution map mapping an initial state x G K to the set S^{x) of solutions to the 
differential inclusion x' G F{x) starting a,t x G K and viable in K. 

The set- valued map S : X C(0, oo; X) is called the solution map associated 
with F. 

We next denote by 

"d^{t,x) := {2;(t)} & := "D^{t,x) 

x(-)es^(x) xec 

the iC-viable reachable maps (or set-valued flow) of x G K and C C K respectively. 
We set d := when viability constraints are absent. 



Definition 31 We associate with the dynamics {F, R) of the impulse differential 
inclusion its substratum F^ := F^ : R+ X K K, that is the set-valued map 
associating with any (t, x) G R+ x K the subset 

F(F^R){t,x) := R{'d^{t,x))riK 

of the elements y G R{c) where c G C := KnR~^{K) through which the solutions 
to the differential inclusion x' G F{x) starting at x and viable in K until they 
reach R~^{K) at time t. 

Knowing the substratum F^ , we introduce 

1. the impulse map 

:= {t > 0 such that F^ j^-^{t,x) yf 0} 

2. and the initialization map : K X 

First, we single out the following property: 

Proposition 32 Let {F, R) be an impulse differential inclusion defined on a 
subset K. Knowing the substratum F^^^^ of {K, F, R), and thus the impulse 
map and the initialization map we can reconstruct a viable run of 

the impulse differential inclusion {F, R) through the following algorithm: Given 
the cadence Tn and the initial state cc„, we take 

( i) the next cadence Tn+i G 

ii) the next reinitialized state x„+i G F^ Xn) C jy{xn), 

Hi) the next motive Xn{-) '■= x(- -h t„), a solution to x' G F(x) satisfying 
^n(O) — ^ G R (^Xji-\-i) 



( 1 ) 
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In other words, in terms of impulse times, given the impulse time tn and the 
initial state Xn, we take 

' i) the next impulse time G 

ii) the next reinitialized state Xn+i G R)i^n+i — tn,Xn) C Ufp ^^{xn) 

< Hi) V t G [tn,tn+i], a solution x{-) to the differential inclusion x' G F(x) 
starting from Xn at time viable in K until it reaches 
at time tn+i- 

( 2 ) 

Proof — Take any run a;(-) associated with a sequence T{x{-)) := {t„} 
of impulse times starting at xq G K and viable in K. Then the sequence x : 
n — >■ x{tn) is a solution of the discrete dynamical system ^-^{tn+i — tn,Xn), 
obviously viable in K. 

Conversely, assume that the substratum F^ is known. The above algo- 
rithm 0 starting at time 0 and state Xq G K provides a run a;(') associated 
with the sequence T{x{-)) := {t„} of impulse times of the impulse differential 
inclusion (F, R) viable in K. □ 

Actually, if we are interested only in the sequence of reinitialized states and 
not necessarily in knowledge of the sequence of impulse times, the knowledge of 
the initialization map is sufficient: 

Proposition 33 A subset K is viable under the impulse differential inclusion 
{F, R) if and only if the domain of the initialization map is equal to K. 

Proof — Assume that K is viable under (F, R) and prove that for every 
X G K, U^p p^{x) 0. Take any xq G K. By definition, there exists a run x(-) 

associated with a sequence T(x(-)) := {t„} of impulse times viable in K. Then 
the sequence x : n ^ x{tn) is a solution of the discrete dynamical system py 
obviously viable in K. 

Conversely, assume that K is viable under the discrete system U(F jiy i.e., 
that for every x G K, U^p-^(x) 0. We shall prove that K is viable under the 

impulse differential inclusion (F,R). Let xq given in K and a solution a; : n — >■ 
Xn G Ufp p^{xn-i) n F to the discrete dynamical system Ufp py By definition 
of the initialization map U^p p -^ , we can associate with 

Xn G Uf^pp^(Xn-l) '■= F^p p^(t, Xn-l) 



some T„_i G such that 

Xn • — Xn(rn—l) G R(^d (^Tn—l,Xn—l)) 

where Xn(-) is a solution to the differential inclusion x' G F(x) starting at time 
0 from x„_i. Setting := tn-\ +r„_i and x{t) := x„(t + t„_i) if t G [tn-iffn], 
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we have checked that x(-) is a run to the impulse differential inclusion (_F, R) 
associated with the sequence {tn}n>o of impulse times tn starting from Xg and 
viable in K. □ 

4 Some Prerequisite from Viability Theory 

Most of the results of viability theory are true whenever we assume that the 
dynamics is Marchaud: 

Definition 41 (Marchaud Map) We shall say that F is a Marchaud map if 

{ i) the graph of F is closed 
ii) the values F{x) of F are convex 
Hi) the growth of F is linear: 3c>0|Va:SX, 

||F(x)|| := sup„g^(^) Ihll < c(||a;|| + l) 

This covers the case of Marchaud control systems where {x^u) i-T f{x,u) is 
continuous, affine with respect to the controls u and with linear growth and 
when P is Marchaud. 

We recall the following version of the important Theorem 3.5.2 of Viability 
Theory, ^ Aubin]: 

Theorem 42 Assume that F : X X is Marchaud. Then the solution map S 
is upper semicompact with nonempty values: This means that whenever Xn € X 
converge to x in X and Xn(-) C S{x„) is a solution to the differential inclusion 
x' G F{x) starting at Xn, there exists a subsequence (again denoted by) Xn(-) 
converging to a solution x{-) G S{x) uniformly on compact intervals. 

Our purpose is to characterize the viability of a subset K under an impulse 
differential inclusion: 

Definition 43 We shall say that a subset K is viable under an impulse differ- 
ential inclusion (F, R) if from any initial state x of K starts at least one run 
viable in K . 

The Viability TheoreirH and its consequences imply the following 

Theorem 44 Let (F, R) be an impulse differential inclusion and K <Z X be a 
closed subset. Assume that F is Marchaud and that R~^{K) is closed. Then the 
following statements are equivalent 

1. K is viable under (F,R), 

2. The K\R ^{K) is locally viable under F, 

^ See for instance Theorems 3.2.4, 3.3.2 and 3.5.2 of P] Aubin]. 

^ The subset K\C denotes the intersection of K and the complement of C, i.e., is the 
set of elements of K which do not belong to C. 
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3. K , F and R are linked through the tangential conditio^ 

Va; G F(a:) n Tjf (a;) ^ 0 

(see |21 Aubin] or im Aubin, Lygeros, Quincampoix, Sastry & Seube] for a 
proof.) 

We shall also need some other prerequisites from Viability Theory: 

Definition 45 Let C G K G X be two subsets, C being regarded as a target, 
K as a constrained set. The subset Capt^(C') of initial states Xq G K such 
that C is reached in finite time before possibly leaving Kby at least one solution 
x{-) G 5(a:o) starting at Xg is called the viable-capture basin of C in K. A subset 
K is a repeller under F if all solutions starting from K leave K in finite time. 
A subset D is locally backward invariant relatively to K if all backward solutions 
starting from D viable in K are actually viable in K . 

We shall use the following characterization of capture basin (see ^ Aubin]): 

Theorem 46 Let us assume that F is Marchaud and that the subsets C G K 
and K are closed. Lf K\C is a repeller (this is the case when K itself is a 
repeller), then the viable- capture basin Capt*^(C) of the target C under S is the 
unique closed subset satisfying C G D G K and 

J i) D\C is locally viable under S 

1 ii) D is locally backward invariant relatively to K ^ 



5 The Graph of the Substratum 

We begin by characterizing the graph of the substratum F^ : 

Theorem 51 Let us assume that F is Marchaud, that C G R is closed and that 
the graph of R : C X is closed. 

Then the substratum F^ : K K is the unique set-valued map with 
closed graph satisfying 

\/xGK, F^j^^{0,x) := R{x)nK 



and, for any T > 0 

1. for any y G F^ j^^{T,x), there exists a solution x(-) to the differential inclu- 
sion x' G F(x) viable in K on [0,T] such that 

vtG[o,r], y^r^p,,^{T-t,x{t)) (4) 

The contingent cone Tl{x) to L G X at x € L is the set of directions v € X 
such that there exist sequences > 0 converging to 0 and Vn converging to v 
satisfying x -\- h„Vn G K for every n (see for instance 0 Aubin & Frankowska]) or 
cni Rockafellar & Wets] for more details). 
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2. for any y G K\F^ , x) , for every solution x{-) to the differential inclu- 

sion x' G F{x) viable in K on [0,T], then 

VfG[0,T], yGK\rlf^^-^{T-t,x{t)) 

As a consequenc^, for any T > 0 and for any y G dxF^ j^-^{T,x), for every 
solution x{-) to the differential inclusion x' G F(x) satisfying then 

VtG[0,r], yGdKF^p^^iT-t,x{t)) 

For proving Theorem we shall first observe that the graph of the substra- 
tum of (K,F,R) is a viable-capture basin and next, deduce the above results 
from the characterization of viable-capture basins. Let us recall that we denoted 
by the graphical restriction of R to K x K defined by 

— f n iC if X G iC 

We observe that C := Dom(i?j^) = iC n R~^{K), that Im(i?j^) and that 
Graph(i?j^) = Graph(i?) C\{K x K). 

Lemma 52 The graph of the substratum F^ of(K, F, R) is the viable- capture 
basin o/ {0} x Graph(i?j^) under the set-valued map {—1} x F x {O}.' 

Graph(T(^_«)) = Gaptf_+"j^;f{„j ({0} x Graph(i?j;^)) 

andy xG C := KnR~\K), F^j^^{0,x) = R{x)nK. 

Proof — Indeed, to say {T,x,y) belongs to the viable-capture basin 

Gapt^_+ ({0} X Graph(i?j^)) 

means that there exists a solution cc(-) G 5(x) and t G [0,T] such that 

ii) Vf G [0,t], {T-t,x{t),y) G Capt^+^j^^f^gj({0} x Graph(i?j^) 

\ ii) (T - t, y, x{t)) G {0} x Graph(i?j^) 

i.e., if and only if f = T and 

(i) ytG [0,T[, x{t) G K 
\ii)y G R{x{T))r^K 

This is equivalent to say that y G F^ x) fl K. 

The relative boundary 8kD to A" of a subset D G K is equal to D fl K\X. 
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Consequently, to say that y belongs to x) means that y G Ji(x)nK. 

□ 

Proof of Theorem 1511 — We observe first that the map {— 1} x P x {0} : 
RxXxX R X X X X is Marchaud and that R+ x K x K is & repeller under 
this map since any solution {T—t, x{t), y) starting at (T, x, y) leaves R+ xKxK 
at time T. Theorem Ei states that the viable-capture basin 

Graph(T(^_^)) = Capt^+^j^^f^gj ({0} x Graph(iij^)) 

is the unique closed subset V d R x K x K containing {0} x Graph(i?j^) 
satisfying 

1. V\({0} X Graph(i?j^)) is locally viable under {—1} x F x {0} 

2. and 

, R,-l- X K X K /•» %\ « 

Capt{_+,j,,^^{oj(V) = V 

This states that whenever (T,x,y) G (R+ x K x K)\V, all solutions to the 

differential inclusion (t',x',y') S {—1} x F{x) x {0} leave (R+ x K x K) 

before possibly reaching the target {0} x Graph(i?j^). 

The first statement means that whenever (T, x, y) belongs to V, there exists 
a solution a;(-) to the differential inclusion x' G F{x) such that (T — t,x{t),y) 
belongs to V until it reaches {0} x Graph(i?j^). This is equivalent to saying that 

VtG[0,r], V & — t,x{t)) 

The second statement means that whenever (T, x, y) does not belong to V, all 
solutions x{-) to the differential inclusion x' G F{x) are such that {T — t,x{t),y) 
do not belong to V whenever (T — t, x{t),y) G R+ x K x K, i.e., whenever x(-) 
is viable in K on the interval [0,T]. This is equivalent to saying that for all 
solutions to x' G F{x) viable in K on the interval [0,T], 

VtG[0,T], yGK\F^pp){T-t,x{t)) 

Let us consider now y G dF^ p<^{T,x) where T > 0. This means that there 
exists a sequence yn G K such that yn G K\F p^{T , x) . Hence (T, does 
not belong to the capture basin of {0} x Graph(i?j^) viable in R+ x K x K. 
Therefore we know that for any solution x(-) G S{x) viable in K on [0,T], 
for any t G [0,T], G K\F^p-^{T — t,x{t)) and, in particular, that G 
K\F^ p-^{0,x{T)) = R{x{T)). Taking any solution x{-) G S{x) satisfying @ 
and the limit when n — >■ -l-oo, we infer that 

VtG[0,r], yGdKF^pp^iT-t.xit)) 

and that 

yGdKR{x(T)) 
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6 Hamilton- Jacobi Characterization of the Substratum 



Before stating the general result characterizing the substratum as a solution to 
a system of first-order partial differential inclusions, let us consider the following 
particular case: 

Proposition 61 Let us assume that f : X X is Lipschitz, r : X X is 
single-valued and continuous, that is continuous, that K is viable under 

(f,r) and differentiable. Then it is the unique solution to the system 

of first- order partial differential equations 



V X € K\C, y j = 1, . . . ,n, 



duj{t, x) 

Wt 






duj(t, x) 
dxi 



fi{x) 



= 0 



or, in a more compact form, 

V.eK\C, + = 0 

satisfying the condition 



y X £ C, u{0,x) = r{x) 



We shall deduce from Theorem El below. Indeed, thanks to the concepts of 
contingent derivative, we shall show that the substratum is the unique 

(set- valued) solution in the “Frankowska sense” to the “Hamilton- Jacobi inclu- 
sion” 



dV(t,x) 

Wt 



dV(t,x) 

dx 



F{x) 



( 5 ) 



satisfying the condition 



Vx e C, V{0,x) = R{x)nK 

We refer to |5I7I Aubin], |0l Aubin & Frankowska] and their references for set- 
valued solutions to systems of Hamilton- Jacobi inclusions. For that purpose, 
we recall that the (graphical contingent) derivative of a set- valued map V : 
R+ X K K may be defined by the relation 



Graph(W(r,x,?/)) := 7Graph(y)(C y) 



Definition 62 We shall say that a set-valued map V : R+ x K K is a 
Frankowska solution to the Hamilton-Jacobi system of first-order partial differen- 
tial inclusions (E)) satisfying the initial condition V{0,x) = i?(x) if its graph is 
closed, if 



V t > 0, y y £ V{t,x), 3 X £ F{x) such that 0 £ DV{t,x,y){—l,v) 
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and if for every v € F{x) 

y t > 0, y y G V{t,x), OG DV{t, x,y){l, —v) 
or 

( i) —vG Tx\k{x) if a; G dK 
ii) —V G Tx{x) if y G dK 

Theorem 63 Let us assume that F is Marchaud, that C := K C\ R~^{K) is 
closed and that the graph of R : C K is closed. 

1. The substratum F^ : K K is the largest set-valued map V : R+ xK 
K with closed graph contained in K x K satisfying 

y t > 0, y G V(t,x), 3v G F(x) such that 0 G DV{t,x,y){—l,v) 

and the condition V (0, x) = R(x) fl K, 

2. If furthermore, F is assumed to be Lipschitz, the substratum F^ ^ ■. K K 
is the unique Frankowska solution V : R+ x K K to the Hamilton- Jacobi 
system of first-order differential inclusions & satisfying the initial condition 
F(0,a;) = R{x). 

Proof — When F is Marchaud, to say that Graph(P)\({0} x 
Graph is locally viable under {—1} x F x {0} means that 

V (t,x,y) G Graph(l/)\({0} x Graph , 

{-1} X F{x) X yf 0 

We observe that (t,x,y) G Graph(P)\({0} x Graph(i?j^)) whenever f > 0 and 
we recall that 

^Graph(y)(^’^>y) = Graph(Hl/(t, x, y)) 
so that the above condition reads 

y t > 0,y y G F^j^^{t,x), 3 V G F{x) such that 0 G DV{t,x,y){—l,v) 
When F is assumed to be Lipschitz, to say that 

Captf_+^j^;f{oj(Graph(G)) = Graph(y) 

means that 

1. for any (t,x,y) G Graph(M)) nlnt(R+ x K x K), 

({1} X -F(x) X {0}) C ?Graph(y)(^’^-y) = Graph(Hl/(t, x, y)) 

This is equivalent to say that for every v G F{x), 

V t > 0, X G Int(iG), y G G(t, x) n Int(AT), 0 G DV{t,x,y)(l,—v) (6) 



The Substratum of Impulse and Hybrid Control Systems 117 



2. and otherwise, for any (t,x,y) G Graph(t/)) fl 9(R+ x K x K), 

({1} ^ ^ {0}) G "^Graphlvi^^’ 2/) 2;, y) 

This means that for every v G F(x), 

( i) 0 G DV{t,x,y){l,—v) lit = 0,y G R{x) 

I ii) 0 G D(t, X, y)(l, —v) or — v G Tx\k > 0,x G dK, y G R{x) 

in) 0 G D{t, X, y)(l, —v) or — v G Tk if t > 0, y G R{x) fl dK 

Indeed, 

(R X X X X)\(R+ xKxK) = 

(R_ xK xK)\J (R+ X {X\K) x K)\J (R+ xKx {X\K)) 

Therefore, condition (1, —v, 0) belongs to the contingent cone to R_ xKxK 
at (0,a:,y) is impossible, condition (1,— 1 !, 0 ) belongs to the contingent cone 
to R^ X (X\K) X K a,t {t,x,y) when x G dK means that —v belongs to 
Tx\k{x) and condition (1,— ?;,0) belongs to the contingent cone to R_ x 
K X (X\K) at (t,x,y) when y G dK means that —v belongs to Tk{x). □ 

For the initialization map, we obtain the following Hamilton- Jacobi inclu- 
sion : 

Theorem 64 Let us assume that F is Marchaud, that C := K C\ R~^{K) is 
closed and that the graph of R : C K is closed. 

1. The initialization map : K K is the largest set-valued map V : 

R+ X K K with closed graph contained in K x K satisfying 

V y GV{x), 3v G F{x) such that 0 G DV{x,y){v) 

2. If furthermore, F is assumed to be Lipschitz, the initialization map Uf^p : 
K K is the unique Frankowska solution V : R+ x K K to the 
Hamilton- Jacobi system of first-order differential inclusions (Ej) satisfying 
the condition x G C, V{x) = R{x). 
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Abstract. Path-dependent impulse differential inclusions, and in partic- 
ular, path-dependent hybrid control systems, are defined by a path- 
dependent differential inclusion (or path-dependent control system, or 
differential inclusion and control systems with memory) and a path- 
dependent reset map. 

In this paper, we characterize the viability property of a closed subset of 
paths under an impulse path-dependent differential inclusion using the 
Viability Theorems for path-dependent differential inclusions. 

Actually, one of the characterizations of the Characterization Theorem 
is valid for any general impulse evolutionary system that we shall defined 
in this paper. 

Keywords: hybrid control, impulse control, path-dependent differential 
inclusion, differential inclusion with memory, functional differential in- 
clusions, viability, run, execution, Kakutani Theorem, contingent cone, 
Marchaud map. 



Introduction 

In this paper, we characterize the viability property of a closed subset of paths 
under an impulse path-dependent differential inclusion using the method of 
[0 Aubin] or 0, Aubin, Lygeros, Quincampoix, Sastry & Seube], the Path- 
Dependent Viability Theorems of 1 1 2f 1 1 4| Haddad]. 

Actually, one of the characterizations of the Characterization Theorem is true 
for any general impulse evolutionary system that we shall define in this paper, 
which is based on recent results of |Sl Aubin]. 

We recall that hybrid control system^ can be embedded in the framework 
of impulse differential inclusions; in the same way, path-dependent hybrid sys- 
tems can be regarded as instances of viable path-dependent impulse differential 
inclusions, and enjoy the same properties. 

^ See for instance among many papers and books [101 Branicky, Borkar & Mitter], 0 
Bensoussan & Menaldi], Matveev & Savkin] and CHI Shaft & Schumacher]. 
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1 Impulse Path-Dependent Differential Inclnsions 

Let X := R" be a finite dimensional vector space. 



1.1 History Spaces 

We denote by 

n{X) := C(-oo,0;X) 
the history (or memory, path) space . 

It is supplied with the compact convergence topology. We denote by 'H\{X) 
the subset of Lipschitz functions with Lipschitz constant A. 

If K c 'H(X), we set 

VtG]-oo,0], K(r) := {</j(r)}^gK 

Observe that Ascoli’s Theorem states that a closed subset K C 'H\{X) is 
eompact if and only if K(0) := {(p(0)},^gK is bounded, since it is closed and 
equicontinuous (by assumption) and pointwise bounded because, for all f/' G K 
and T < 0, 



mr)\\ < Mr) - m\\ + IIV'(0)|| < A|r| + ||K(0)|| 

Our study invlves a constrained subset K C T~L{X) made of paths or histories 
and of a target C C K. 

A first example of constrained subset of paths or histories and targets asso- 
ciated with subsets C d K CL X oi the vector space are given by 

C := {ipC n{X) I (^(0) G C } C K := { G n{X) I V3(0) G K } 

where the constraints bear only on the present. 

Another class is given by Volterra sets defined through a “kernel” k :] — oo, 0] x 
X 1 —^ Y and a set- valued map M : Y X hy 

K := G ?^(A) I <^(0) G m(^J° fc(-s,(^(s))d/x(s)^| 

where the constraints involve cumulated consequences of the history. 

In the discrete case, 



K := 



G H{X) I <p(0) G 




K-hvU^j)) 



j = -oo 



involves discrete cumulated consequences of the history (delays). 

Associated targets can be asociated with set- valued maps P C M in the same 
fashion. 
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1.2 Histories or Paths 



We associate with any continuous function x(-) € C(—oo, +oo; X) its history ( or 
path^ T{t)x up to time t defined by: 

VtG] — oo,0], T{t)x{T) := x{t + T) 

Then T{t) maps C(— oo, -boo; X) to 'H{X) and satisfies the semi-group property 

T(t + s)x = T{s)T{t)x 

We then observe that for any function x(-) S C(— oo, -boo; X), we have x{t) = 

(T(t)x)(0). 

In this continuous framework, we define the constraints of the history of the 
evolution through a closed subset K C 'H(X). Viable evolutions x{-) with 
memory are the ones that satisfy 

Vt>0, T{t)x e K (1) 



and an evolution x{-) reach a target C at time s if T(s)x € C. 

For instance, an evolution is viable in K :={(/? G T~L{X) \ i^(0) G K} if and 
only if for every t g] — oo, 0], x{t) G K. If 



K := G "H(X) I <p(0) G (/7(s))d/x(s) 

then a:(-) is viable in K if and only if 

V t > 0, x{t) € m( f k(t — s,ip{s))dp,{s) 

\J — OO 



1.3 Path-Dependent Differential Inclusions 

Let us consider a set- valued map F : "H(V) H> V governing the continuous 
evolution of the state x(t) through the path-dependent differential inclusion 

for almost all t > 0, x'{t) G F{T(t)x) 

starting at a given ip G TL{X) in the sense that 

r(0)a; = p 

i.e., for every r g] — oo,0], x{t) = p{t). 

We denote by TZp ■ 'H{X) C(0, oo; X) the map associating with any initial 

path p G 'H{X) the set TZp{p) of solutions t >->• x{t) to the path-dependent 
differential inclusion x'{t) G F{T{t)x) starting at the initial path p in the sense 
that ^(O)^ = p. 

often denoted by xt := T(t)x 



2 
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Actually, we shall need the properties of the associated set- valued map Sp ■ 
T~L{X) ^ C(0, oo; 'H(X)) defined by 

Spif) ■■= {t (p{t) := T{t)x}^(^.)eTZF{v) 



Definition 11 We shall say that this set-valued map Sp ■ T~L{X) C(0,-|-oo ; 
TL{X)) is the solution map of F. 

The solution map Sp has the advantage of mapping the set T~L{X) into time- 
dependent functions t ‘p{f) := T{t)x that belong to the same histoy space 
TL{X), even though the traditional view is to call a solution a function 1 1 — x{t), 
taking its values in X. 

This choice of Sp instead of TZp is justified by its following properties: 



1. the translation property: Let (/?(•) G S((p). Then for all s > 0, the function 'if(-) 
defined by 'ip(t) := (p(t s') is the history ^/>(-) := T(-)y G S(T(s')x) of the 
solution y(-) to the path-dependent differential inclusion starting at T{s)x, 

2. the concatenation property: Let (p{-) G Sp{(f) be the history of a solution to 
the path-dependent differential inclusion starting at the path (p and s > 0. 
Then for every history ip{-) G Sp{T{s)x) of a solution y{-) to the path- 
dependent differential inclusion starting at the initial path T(s)x, the func- 
tion ^(-) defined by 

’ \if{t — s) := T{t — s)y if t > s 
is the history of the solution z(-) defined by 



z{t) 



f x(t) if t G [0, s] 
y(t — T) if t > s 



to the path-dependent differential inclusion starting at the initial path p, 
and thus, belongs to Sp{ip). 



These two properties to which we add the upper compactness of the solution 
map are enough to obtain relevant (and interesting) properties of path-dependent 
impulse differential inclusions, common to other dynamical systems. 



1.4 Runs of Impulse Path-Dependent Differential Inclusions 

We now introduce a constrained functional set K C TL{X), a functional target 
C C K and a path-dependent reset map i? : C ^ K with nonempty values 
R{(p). 

The pair (T, R) governs the evolution of impulse systems in the following 



sense. 




Path-Dependent Impulse and Hybrid Systems 123 



Definition 12 Let us consider a finite dimensional vector space X, the space 
'H{X) of histories, a subset K C 'H{X), a target C C K, a set-valued map 
F : 'H(X) X and a set-valued map R : C ^ Ti with nonempty values, regarded 
as a path-dependent reset map. We regard the pair (F,R) as the dynamics of a 
path-dependent impulse differential inclusion. 

A run of the path- dependent impulse differential inclusion (F,R) is defined 
by 

1. a finite or infinite sequence r(x(-)) := {r„}„ of nonnegative cadences t„ G 
[0,oo[, 

2. a sequence of reinitialized paths ipn G R{X), 

3. a sequence of motives ipn{') '■= T(-)x„ G SF{ipn) where ipn G R{X) is the 
history of a solution Xn(-) to the path- dependent differential inclusion x'{t) G 
F{T(f)x) starting at the initial path ipn o,nd satisfying the end-point condition 
T{Tn)Xn G R~^{ifn+l) 

by 

( i) defining the sequence of impulse times tn+i ■= tn -\- Tn, ,r,\ 

ii') V f G [tn, fn-t-l [j ^(f) ■ — ^nft tjf) 

If the sequence of cadences is finit^ and stops at tn, we agree that the Nth 
motive is defined on [0,-|-oo[, i.e., that we take tn+i = -l-oo. 

We say that a run x{-) is viable in K if for any t > 0, T(t)x G K and that 
K is locally viable under {F, R) if for any ip gK., there exists at least one run 
of the impulse path- dependent differential inclusion viable on a nonempty time 
interval and (globally) viable if it is viable on [0, -|-oo[. 

At this stage, a run x(-) can just be a (discrete) sequence of paths Pn+i G 
R{(fin) at the initial time (case when for all n > 0, the cadences r„ = 0), or just 
a (continuous) solution a;(-) to the path-dependent differential inclusion x'(t) G 
F{T{f)x) (case when t\ = +oo), or an hybrid of these two path-dependent 
modes, the discrete and the continuous. 

Path-dependent hybrid systems can be regarded as instances of viable path- 
dependent impulse differential inclusions as in the case of usual hybrid systems: 
we refer to 0 Aubin] or 0 Aubin, Lygeros, Quincampoix, Sastry & Seube] for 
more details on this topic. 

2 Statement of the Impulse Path-Dependent Viability 
Theorem 

2.1 Marchaud Maps 

The Viability Theorems hold true whenever we assume that the dynamics gov- 
erning the path-dependent evolution is Marchaud: 

^ We shall see that we can eliminate this situation by assuming that R{C) n 
ViabF(K) = 0. 
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Definition 21 (Marchaud Map) We shall say that F : 'H(X) X is a Mar- 
chaud map if 



( i) F is upper semicontinuous 
ii) the values F{ip) of F are convex 
Hi) the growth of F is linear:3 c > 0 | V G ’H(AT), 
||FM||:=sup„,^(^)||^|| < c(||^(0)|| + l) 

This covers the case of Marchaud control systems where {(p, u) i— >■ f{ip, u) is 
continuous, affine with respect to the controls u and with linear growth and 
when P : 'H(X) Z is Marchaud. 

We recall the following version of the important Haddad Theorem 12.4.1 of 
IP Aubin]: 

Theorem 22 Assume that F : 'H(X) X is Marchaud. Then its solution 
map Sp is upper semicompact with nonempty values: This means that whenever 
Pn G T~L{^) converge uniformly on compact intervals to ip in TL{X) and any 
history Pn{') '■= T{-)xn G SpiPn) associated to a solution Xn(-) to the path- 
dependent differential inclusion x'(t) G F(T(t)x) starting at pn, there exists a 
subsequence (again denoted by) Pn(') converging uniformly on compact intervals 
to the history p{-) := T{-)x of a solution x{-) to the path- dependent differential 
inclusion starting at p. 



2.2 Contingent Directions 

In the case of path-dependent impulse differential inclusions, we shall charac- 
terize the viability of a functional constrained set K in terms of contingent 
directions to a K C ?t(A) be a subset of histories at a path p gTL. Let 

A{X) := {a;(-) G C(0, -boo; A) such that cc(0) = 0} 

denote the “future space” . We embed the state space X into A{X) by identifying 
a vector x with the function ifx{t) := tx. The image of the ball XBx of radius A 
under this embedding is contained in A\{X) of A-Lipschitz functions. 



Definition 23 Let h > 0 be given. The ^-concatenation ('or concatenation when 
there are no ambiguities) p'O^ip is the bilinear form from TL{X) xA{X) i— TL{X) 
defined by 



{pOhiffir) 



f p(t -h h) if r g] — oo, —h] 

\ p(0) -h tp{T -b /i) if r G [—h, 0] 



As an example, the concatenation of p G TL{X) and a; G A is defined by 



{pOhx){T) = 



p(t -\- h) if r g] — 00, h] 

p(0) -b (r -b h)x if r G [— /i, 0] 
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We attach to a function ip G 'H(X) and a sequence vi,...,Vn the Euler 
concatenation defined by 



‘p0hVi0hV20h ■ • • OhVn 



which is piecewise linear function on the interval [—nh, 0]. 

Lemma 24 For any Lipschitz constant A > 0, the h- concatenation maps 
HxiX) X Ax{X) to Hx{X) 



Definition 25 We denote by V-K^ip) the set of vectors v G X such that there 
exist a sequence > 0 converging to 0 and a sequence Vn G X converging to v 
satisfying 

V n > 0, ipOh^Vn G K 

2.3 The Impulse Path-Dependent Viability Theorem 

Theorem 26 Let (f, R) be a path- dependent impulse differential inclusion and 
K C R-iX) be a closed subset. Assume that F is Marchaud and that C C K is 
closed. Then the following statements are equivalent 

1. K is viable under (F,R), 

2. The subset K\C is locally viable under the path- dependent differential inclu- 
sion governed by F , 

3. K, C, F and R are linked through the tangential condition 

V(^eK\c, 

Actually, both iimulse differential inclusions and path-dependent impulse 
differential inclusion^ share the same properties at a higher abstraction level, 
the level of impulse evolutionary systems we are about to define. It is at this 
level that the two first statements are equivalent. 

The equivalence between the second and third statement is specific, and 
provided in our case by the Path-Dependent Viability Theorems of miiii 
Haddad] . 



2.4 Examples 

Take any path-dependent differential inclusion x'{f) G F(T(t)x) associated with 
a Marchaud right-hand side F. 

We refer to Chapter 12 of ^ Aubin] for examples of tangential conditions 
when the constrained set K and the constrained targets C C K are Volterra 
sets defined by kernels, by lack of space. Consider only the simple case when 

as well as parabolic (or reaction-diffusion type) partial differential inclusions and 
mutational equations governing the evolution of subset. 
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the constrained subset of paths or histories and the target are associated with 
subsets C d K <Z X oi the vector space X by 

C := { (^ G %{X) I (^(0) G C } C K := { V? G %{X) \ ^^(O) G K } 

Assume that the reset map R is associated with an usual reset map Rq '■ C 
K, where C C K hy the formula 

V V? G "H(A), Vr g] - oo,0], (i?((p))(r) := i?o(v5(T)) 

In this case, a run of the path-dependent impulse differential inclusion (F, Rq) 
is defined by 

1. a finite or infinite sequence r(a;(-)) := {r„}„ of nonnegative cadences t„ G 
[0,oo[, 

2. a sequence of reinitialized states a;„ G K, 

3. a sequence of motives x„(-) that are solutions to the path-dependent differen- 
tial inclusion x'{t) G F{T{t)x) starting at the initial state and satisfying 
the end-point condition Xn+i G R{xn{Tn)) 

by 

J i) defining the sequence of impulse times := -I- Tn, /o', 

\ it) V t G x{t) := Xn(t-tn) 

Theorem 27 Let F be a path-dependent Marchaud set-valued map, K and C C 
K and Ro : C ^ K be a reset map. Then K is viable under {F, Rq) if and only 
if the tangential condition 

y (p G 'H(A') such that (^(0) G K\C, F{ip) nT/f(v3(0)) yf 0 



2.5 Path-Dependent Hybrid Systems 

Definition 28 An path-dependent hybrid differential inclusion {K, F, Rq) is de- 
fined by 

1. a finite dimensional vector space E of states e called locations, 

2. a set-valued map K : E X associating with any location e a (possibly 
empty) subset K{e) d X and a set-valued map C : E X associating with 
any location e a (possibly empty) subset C{e) d K{e), 

3. a set-valued map E : E x 'H(X) X with which we associate the path- 
dependent differential inclusion x'(t) G E(e,T(t)x), 

4- a set-valued map (reset map) Rq : Graph(C) Graph(AT). 

A run of the path- dependent hybrid system is defined by 

1. a finite or infinite sequence r(e,a;(-)) := {t„}„ of nonnegative cadences r„ G 
[0,oo[, 

2. a sequence of locations e„ and of reinitialized states x„ G K{en), 
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3. a sequence of motives a;„(-) that are solutions to the path- dependent dif- 
ferential inclusion x'{t) € F{en,T{t)x) starting at the initial state Xn and 
satisfying the end-point condition Xn+i G i?o(e„, x„(r„)) 

by (0). 

We observe right away that a map (e,a;(-)) is a run of the hybrid differential 
inclusions if and only if (e(-), a;(-)) is a run of 

r t) e'{t) = 0 

G F{e{t),T{t)x) 

"viable” in Graph(itT) until it reaches the graph of the map C. Indeed the loca- 
tions remain constant in the intervals [t„, fn+i[ since their velocities are equal to 

0 . 

Since the existence of solutions to path-dependent hybrid differential inclu- 
sions amounts to the viability of the graph of the set-valued map K under an 
associated auxiliary path-dependent impulse differential inclusion, we obtain a 
necessary and condition for the existence of solutions to hybrid differential in- 
clusions thanks to Theorem For that purpose, we need the definition of the 
contingent derivative DK{e,x) ■. E X of a set- valued map K : E X at a 
point (e,x) of its graph: It can be defined by 

Graph(i:)iF(e,a:)) := Tqj.^^Yi(k)(^,x) 



Theorem 29 Let (K,F,Rq) be a path- dependent hybrid differential inclusion. 
Assume that F is Marchaud. Then the path- dependent hybrid differential inclu- 
sion has a solution for every initial state if and only if 

y e € E,\/ (fi G TL{X) such that 

^{0) G Kie)\K{e)\C{e), F(e,ip) D DK{e,m){0) ^ 9 

3 Impulse Evolutionary Systems 

Therefore, it costs nothing to prove the equivalence between the two first state- 
ments in the general case of impulse evolutionary systems: 

3.1 Impulse Evolutionary Systems 

Definition 31 An evolutionary system is a set-valued map S : X C{0, oo; X) 
satisfying 

1. the translation property; Let xf) G S(x). Then for all T > 0, the function 
yf) defined by y{f) := x{t-\-T) is a solution y{-) G S(x(T)) starting at x{T), 
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2. the concatenation property; Let x{-) G S(x) and T > 0. Then for every 
y{-) G S{x{T)), the function z{-) defined by 



z{t) 



(x{t) iftG[0,T] 
\y{t-T) ift>T 



belongs to S{x). 

We can define impulse evolutionary systems in the following way: 

Definition 32 Let K <Z X , C G K be two nonempty subsets and R : C K 
a set-valued majl with nonempty values, regarded as a reset map, and S : X 
C(0,oo; AT) be an evolutionary system. Then the pair (S,R) governs a run x(-) 
of an impulse evolutionary system defined by 

1. a finite or infinite sequence t(x(-)) := {r„}„ of nonnegative cadences G 
[0,+oo[, 

2. a sequence of reinitialized states Xn, 

3. a sequence of motives Xn(-) G S(xn) satisfying the end-point condition 

G R (Xyi+i) 

by 

J i) defining the sequence of impulse times tn+i '■= tn + Tn, 

^ H) V f G [tn; tn+1 [; ^(t) ■ — Xjii^t tn) 

If the sequence of cadences is finit^ and stops at tn, we agree that the Nth 
motive Xjy{-) G S{xn) is taken on [0,+oo[, i.e., and we agree to setrjv+i = +oo. 

We say that a run a;(-) is viable in K if for any t > 0, x{t) G K and that 
a closed subset K is viable under an impulse evolutionary system {S,R) if from 
any x G K starts at least one run viable in K . 

In order to characterize the viability of K under an evolutionary system, we 
also need the following definitions: 

Definition 33 Let S \ X C(0,+oo; Al) be a set-valued evolutionary system 
and K G X be a subset regarded as a constrained set. 

The subset K is said locally viable under S if from any initial state x G K 
starts at least one solution viable in K on a nonempty interval and viable if this 
solution is viable on [0,+oo[. 

The viability kernel Viab(AT) is the subset of initial states xq € K such that 
one solution x{-) G 5(a;o) starting at Xg is viable in K for all t > 0. A subset K 
is a repeller under S if its viability kernel is empty. 

® When R : X X is defined on X, we associate with it its “graphical restriction” 
to K X K (again denoted by) R where C ~ K C\ R~^{K) and R{x) is replaced by 
R{x)nK. 

® We shall see that we can eliminate this situation by assuming that R{C) n 
XiahsiK) = 0 . 
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3.2 Characterization of Impulse Evolutionary Systems Viability 

We can adapt the proof of Pl Aubin] or Pl Aubin, Lygeros, Quincampoix, Sastry 
& Seube] for characterizing the viability of a closed subset under an impulse 
evolutionary system: 

Theorem 34 Let (5, R) be an impulse evolutionary system and K G X and 
C G K he elosed subsets. Assume that S is upper semieompaet. Then the follow- 
ing statements are equivalent 

1. K is viable under (S,R), 

2. The subset K\C is locally viable under S, 

3.3 Prerequisites of Viability Theory 

For proving this characterization theorem, we need some results of viability 
theory. 

Definition 35 Let K G X he a subset. The functional : C(0,oo;A) i— >■ 
R+ U {-boo} associating with x{-) its exit time defined by 

tk{x{-)) := inf (t G [0,oo[ I x(t) ^ a:} := ti7x\k{x{-)) 

is called the exit functional. 

Let C G K be a target. We introduce the (constrained) hitting functional 
'^(K,c) defined by 

W(^K,c)i.x{-)) := infjt > 0 | x{t) G C & Vs G [0, t], a;(s) G AT } 

associating with x{-) its hitting time, introduced in m Cardaliaguet, Quincam- 
poix & Saint-Pierre]). 

Consider an evolutionary system S : X C(0, -boo; X). Let C G K and K 
be two subsets. 

The function \ K ^ R+ U {+oo| defined by 

Xk(^) ■= sup TKix{-)) 

x{-)^S(x) 

is called the upper exit function. 

The function : K^R+ U l+oo} defined by 

is called the lower constrained hitting function. 



We shall need the following 
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Theorem 36 Let S : X C(0,+oo;X) be a striet upper semieompaet map 
and C and K be two closed subsets such that C <Z K . Then the hitting function 
C) semicontinuous and the exit function is upper semicontinu- 

ous. Furthermore, for any x € Dom(zu^j^ there exists at least one solution 
x^{-) € S{x) which hits C as soon as possible before possibly leaving K 

and for any x G Dom(r^), there exists at least one solution G S{x) which 
remains viable in K as long as possible: 

Tk(^) = tk{xH')) 

(See 13 Aubin] for a proof and more details on evolutionary systems). 

3.4 Proof of the Characterization of Impulse Evolutionary Systems 
Viability 

Indeed, if K is viable under (S,R), then from any xq G K\C starts at least a 
solution x{-) G S{x) viable in K 

1. either forever if xq belongs to the viability kernel Viab(A') of K 

2. or until it reaches at some time > 0 a state x{~ti) in C. 

This shows that K\C is locally viable. 

Conversely, let us assume that K\C is locally viable and take an initial state 
Xq G K. If Xq belongs to C, we may take tq = 0 and x\ G R{xef). Consider now 
the case when xq G K\C. 

If Xq belongs to Viab(AT), then at least one solution starting from xq is viable 
in K, and thus, defines a run viable in K: We may take the cadence tq = +oo 
and for motive a solution Xo(-) G S(xq)- 

If xo does not belong to Viab(AT), all solutions leave K in finite time before 
(possibly) reaching the viability kernel. It is then enough to prove that at least 
one of them reaches C before leaving K. This is the case of a solution x^(-) G S(x) 
which maximizes tk{x{-)), i.e., which satisfies 

rj^(x) := sup tk{x{-)) = tk(x^(x)) 

x{-)^S{x) 

leaves Ar\(Viab(AT) U C) through C. This solution exists by Theorem EEl since 
K is closed and S is upper semicompact. Next, we claim that x'^ := x**(r^(x)) G 
K\Viab{K). Otherwise, if x^ would belong to the viability kernel, it could be 
concatenated with a solution viable in K for ever, so that the initial state xq 
would belong the viability kernel, which is not the case. 

Furthermore, x^ belongs to C. If not, x^ would belong to K\C which is 
assumed to be locally viable. Then one could associate with x® G K\(Viah{K) U 
C) a solution y(-) G 5(x*) and T > 0 such that ?/(t) G Ar\(Viab(AT) U C) for all 
r G [0, Tj. Concatenating this solution to x**(-), we obtain a solution viable in K 
on an interval [0,r^(x) + T], which contradicts the definition of x**(-). 

Therefore x^ belongs to Ff fl C so that there exists x{ G K (1 R{x^). □ 
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4 The Path-Dependent Viability Theorem 

Since the solution map of a Marchaud map F : 'H(X) ^ X is upper semicompact 
by Theorem E21 the equivalence between the first and second statements of 
Theorem El holds true. 

The equivalence between the second and the third statement follow from the 
following Haddad’s Path-Dependent Viability Theorem: 

Theorem 41 Assume that F is Marchaud and take A > 0. The two following 
statements hold true: 

1. IfK C 'H\{X) is closed, then K is (globally) viable under F if and only if 

V(^eK, ^0 

2. If C cJA is closed, then K\C is locally viable under F if and only if 

V:^eK\C, F{g:)nVK{<f)^0 
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Abstract. In this paper, we consider the problem of stabilizing the kine- 
matic model of a car to a general path in the plane, snbject to very mild 
restrictions. The car model, although rather simplified, contains some of 
the most relevant limitations that make application of existing results 
in the literature impossible: namely, the car can only move forward, and 
turn with a bounded steering radius; also, only limited sensory informa- 
tion is available. 

The approach we follow to stabilization is to adapt to the present gen- 
eral case an optimal synthesis approach successfully applied in our pre- 
vious work to tracking rectilinear paths. Due to both the nature of the 
problem, and the solution technique used, the analysis of the controlled 
system involves a rather complex switching logic. Hybrid formalism and 
verihcation techniques prove extremely useful in this context to formally 
proof stability of the resulting system, and are described in detail in the 
paper. 



1 Introduction 

In this paper we consider the design of a control law for path tracking by a 
so-called Dubins’ model of a car. Dubins’ cars are kinematic models of wheeled 
(nonholonomic) vehicles that move only forward in a plane, and possess a lower- 
bounded turning radius. The model is relevant to the kinematics of road vehicles 
as well as aircraft cruising at constant altitude, or sea vessels. 

Although the design of control techniques for nonholonomic vehicles has been 
the subject of extensive research recently (see e.g. mm), the additional con- 
straint that the steering radius of the vehicle is lower bounded has not been 
explicitly considered. However, such a restriction appears to be crucial in mak- 
ing a kinematic model of a car relevant to real-world vehicles encountered in 

* The work has been conducted with partial support of PARADES, a Cadence, 
Magneti-Marelli and ST-microelectronics E.E.I.G, by CNR PF-MADESSII SP3.1.2. 
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most applications. Another important assumption often used in the literature is 
that the full state of the system is available for measurement, and that the path 
to be tracked is entirely known in advance. Instead, we consider in this paper 
the more realistic and less demanding case that the vehicle can only measure 
its current distance and heading angle error with respect to the closest point 
on the reference path in the plane, where only the sign of the path curvature is 
detected. 

The approach we follow to stabilization of Dubins’ cars is to adapt to the 
present general case an optimal synthesis approach successfully applied in our 
previous work to tracking rectilinear paths m- Due to both the nature of the 
problem, the type of sensors, and the solution technique used, the analysis of the 
controlled system involves a rather complex switching logic. Hybrid formalism 
(see mm) and verification techniques (see |iSI7ll| l prove extremely useful in 
this context to formally proof stability of the resulting system, and are described 
in detail in the paper, which is organized as follows. 

In Section 0 a hybrid automaton that describes the motion of the vehicle 
with respect to the path is introduced, while in Section 0the path-tracking con- 
troller is developed. Such controller, described in detail in Section llOl is obtained 
by considering a local approximation of the desired path with the tangent line, 
and by using a feedback controller designed for stabilization on straight paths 
(reported in Section 13.11 . The advantages of the novel hybrid path-tracking 
formalization are exploited in Section 0 where the stability properties of the 
proposed controller are investigated. By a reachability analysis in the continu- 
ous state space, a finite state abstract representation of the hybrid closed-loop 
automaton is obtained. Though this representation is not a bisimulation, but 
rather a simulation, of the hybrid automaton (| 5 |), it suffices to prove the sta- 
bility properties of the proposed control. It is shown that the proposed hybrid 
feedback controller achieves stabilization of the Dubins’ car on a generic reference 
path and sufficient conditions for global attractivity are derived. 

2 Hybrid Path Tracking Modeling Using Switching 
Prenet’s Frames 

We consider the kinematic model of a car moving forward on a plane, which 
was introduced by Dubins in A configuration of the vehicle is defined by 

an ordered pair {M{x,y),9) G x S^, where {x,y) are the coordinates of a 

reference point M in the plane and 9 is the angle made by the direction of the 
car with respect to the a;-axis. The kinematics of the car are described by 




( 1 ) 



where V is the constant forward velocity, tu the is turning speed and the input 
constraint models a lower bound i? > 0 on the turning radius of the Dubins’ car. 
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Fig. 1. Reference path and transformed coordinates. 



The problem we are concerned with is that of steering the vehicle to a given 
feasible path F, defined in the arclength parametrization by 

r = {(i,y) G I (x,y) = g{/3) for /3 G R} , (2) 

with the following conditions: 

A) g{-) is a class mapping from R to R^ and the orientation of T is that 
induced by increasing. /3; 

B) Let k(/ 3) denote the extension by continuity from the lefi0of the curvature 
of r, expressed as a function of the arclength /3. There exists a positive real 
Rp such that the normalized curvature k{s) = Rk{s) satisfies 

\k{P)\ = R\4P)\<^ = C<l. (3) 

C) Considering the open neighborhood of the path 

Tr = {x G R" : 3/3 e R, ||x - g{(3)\\ < Rp} C R^ (4) 

for all X C 7r there exists a unique nearest point on F. 

In order to describe the motion of the vehicle with respect to the reference path T 
a mobile Frenet’s frame associated to the curve F is considered. Given a vehicle 
position M{x,y) G Tp, the Frenet’s frame is defined by the tangent, 

the principal normal and the binormal axes of the curve at the point j/(/3)) 

of F, located at the minimum distanc^ from M{x,y) (see Figure As the 
vehicle moves with velocity V, the Frenet’s frame St{s) follows its motion so 

^ By definition, k(/ 3) = limg_>^- k(s), at points (x{f3), y{P)) where the curvature of F 
is not defined. 

Note that, by A), B) and C) the Frenet’s frame is well-defined along F. 



2 
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downwards upwards 

k{s) < 0 




Fig. 2. Hybrid automaton PTHA^ modeling the car in the transformed state space. 



as to keep it on the principal normal axis. The arclength abscissa s locates the 
current Frenet’s frame. The tangent and the principal normal axes of St(s) 
remain within the plane containing the curve, while the binormal axis points 
either upwards, if the local curvature of T is counterclockwise (i.e. k(s) > 0), 
or downwards, if the local curvature is clockwise (i.e. k(s) < 0). Introduce the 
transformated coordinates (s,y,0), where: 

— abscissa s defines the position of the Frenet’s frame along the curve; 

— y denotes the position of the car along the principal normal of St{s) (lateral 
distance) normalized with respect to the minimum turning radius R; 

— 9 denotes its orientation with respect to the tangent axis of St{s) (heading 
angle error), with sign taken according to the local direction of the binormal 
axis (see Figure QJ). 

It can be noticed that this coordinate system is similar to the one used by 
Samson jOj, except for the switchings of the Frenet’s frame. In fact, a change 
of curvature along the path produces a jump of the variables y and 9 to the 
symmetric point with respect to the origin in the (y,0)-plane. The reason for 
introducing such discontinuity in the model is related to the different behaviors 
that a vehicle with bounded curvature has when it approaches a reference path. 
Indeed, the approach is apparently easier if the vehicle and the center of cur- 
vature of the path lie on the opposite sides of the curv^. This formulation will 
turn out to be useful in the verification of the proposed path tracking controller. 

The motion of the car in the transformed state (s, y, 0)^ can be described by 
using the formalism of hybrid automata (see m)- The discrete nature of the 
model arises from the fact that the Frenet’s frame St{s) changes its orientation 
during the motion, depending on the sign of the curvature k(s). The discrete 
state, referred to as bin , models the two possible orientations of the binormal 
axis of ST{s(t)) at time t and assumes either the value upwards or the value 

^ For instance, if the vehicle is required to approach a circle with curvature 1/R, then 
it can approach it only from outside. 
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downwards upwards 

ar = switch 




Fig. 3. Hybrid automaton PTHA2 of the vehicle in the reduced state space. 



downwards. Its initial value is: upwards, if k(s(0)) > 0; downwards, if k(s(0)) < 0; 
and any of those, otherwise. The dynamics the continuous states are subject 
to are obtained by geometric arguments. The complete Path-Tracking Hybrid 
Automaton, referred to as PTHAg , is depicted in Figure El 

The specification for the design of a path tracking controller for the Dubins’ 
car can be formulated using the hybrid automaton PTHAg , which captures the 
different behaviors of the bounded-curvature vehicle in approaching the path. 
For such hybrid model, the problem reduces to that one of steering {y, 6) to 
( 0 , 0 ). 

Assuming that only the sign of k{s) is available but not its amplitude, a 
reduced hybrid automaton can be considered for the path tracking problem. 
The local curvature |k(s)| is replaced by an unknown input disturbance d(t) the 
path tracking controller has to be robust to. By 0 , disturbance d{t) satisfies 

0 < d{t) <C <1. (5) 

The path tracking problem is described in the reduced continuous state space 
{y,9). Curvature sign switching conditions k{s) > 0 and k{s) < 0 are modeled 
by a discrete uncontrollable input a^, assuming either the value switch (when a 
change of curvature sign occurs) or the silent move e (otherwise). The reduced 
hybrid automaton, referred to as PTHA 2 , is reported in Figure 0 

In this case the path tracking problem is formulated as follows: 



Problem 1. Let P as in ( 0 ) be a feasible reference path. Given the hybrid au- 
tomaton PTHA 2 , find a feedback control law u)(hin, (y, 0)) satisfying curvature 
constraint (0) such that, from any initial state {bino,{yo,9o)) the trajectory 
(y(f), 9{t)) converges to the origin under the action of any unknown disturbance 
d{t), bounded as in 0, and any sequence of uncontrollable events ar- 
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Table 1. Partition of domain used to define the shortest path synthesis. 




3 Hybrid Path— Tracking Feedback Controller 



3.1 Optimal Feedback Control for Line Tracking 

In dH, the problem of driving the Dubins’ car to a straight path has been 
considered. An optimal feedback control that minimizes the length travelled by 
the vehicle to reach the specified path was deviced. Define aN{y,0) = y + I + 
cos{0) and <jp{y,0) = y — 1 — cos{0). The optimal feedback control presented 
in El is defined inside the region 



V 



iv,e) 



' (JN{y,0) <0 A 0 G [tt, Itt) V 
o-p{y,0) < 0 A 0 e (f ,7 t) V 
< 0G[-f,f] V 
CTAr(y,6») > 0 A 0 G [-7T, -|) V 
. crp{y, 0) > 0 A 0 G (-§7r, -tt) 



( 6 ) 



in the state space {y, 0), which, modulo 2tt angles on 0, corresponds to the whole 
space (see Figure 0. The optimal controller is described by three modes, 
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• go-straight, where a; = 0 

• turn-right, where w = (7) 

• turnJeft, where uj = + ^ , 

which are chosen as follows 

[go-straight, if {y, 6) S A [turn_right, if (y, 0) G f2~] A [turnJeft, if (y, 9) G 17'*"] 

( 8 ) 

where the partition 17° U 17 U 17+ of domain is defined as in Table [D 

In Figure 0the boundaries between the subsets of the partition 17° U 17“ U 17+ 
are represented by dotted lines, and the direction of motion, when the reference 
path is a straight line i.e. d = 0, is represented by directed curves. 

3.2 Feedback Tracking Control for Generic Path 

In this section a hybrid feedback controller that solves Problem Q] is derived 
from the one reported in the previous section. The hybrid model of the vehicle 
PTHA 2 is characterized by the two modes: upwards and downwards. In mode 
downwards input lo appears with opposite sign with respect to mode upwards. 
Since the controller modes in ® has been set assuming an upwards binormal 
axis then, the controller modes turn-right and turnJeft have to be switched when 
the vehicle is in mode downwards. Hence, for a generic feasible path T, the full- 
state feedback controller is defined in {upwards, downwards} x by setting 

the controller modes as follows 

• go-straight, if {bin, (y, 9)) G {upwards, downwards} x 17° 

• turn-right, if {bin, (y, 9)) G {upwards x 17“) V {bin, (y, 9)) G {downwards x 17+) 

• turn-left, if {bin, (y, 9)) G {upwards x 17+) V {bin, (y, 9)) G {downwards x 17“) 

( 9 ) 

where 17°, 17 and 17+ are as in Table Q] The closed-loop hybrid automaton 
CLHA obtained by applying the feedback 0,0 to the vehicle hybrid automaton 
PTHA 2 is depicted in Figure ^ According to 0 and 0 , CLHA has a discrete 
state mode that assumes values in the set O = {zero, negative, positive}, as follows 

• mode = zero if (y, 9) G 17° 

• mode = negative if (y, 9) G 17“ (10) 

• mode = positive if (y, 9) G 17+ . 

The initial state {modeo, {yo,9o)) of the hybrid automaton CLHA has to sat- 
isfy (ITUl) . 

The coordinate transformation {x,y,9) — >■ {s,y,9) becomes singular when 
the vehicle lies on the center of the local osculating circle to the path P. That is 

if, at some time t, y{t) |k(s(t))| = 1, or equivalently y{t) d{i) = 1. For any initial 

configuration {M{xo,yo),So)^ with M{xo,yo) G 7r as in 0, the corresponding 
state {yo,9o) satisfies yo < C~^. Further, since by 0 d < (7, then yod < 1 at 
the given initial condition. However, to ensure that 



y d < 1 i.e. 1 — y d > 0 



( 11 ) 
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Fig. 4. Hybrid model of the closed-loop system CLHA . 

will hold along all the trajectories of CLHA , we need to further restrict the 
admissible initial vehicle configurations, in terms of its initial orientation 9q. 

Proposition 1. Let the continuous disturbance d be bounded to belong to the 
interval [0, C], with 

C < 0.5 . (12) 

Then, m is satisfied along all trajectories of CLHA provided that the initial 
configuration (modco, (yo,0o)) is such that 

(yo,0o) e = |(y,0) e I \y\ < C~^ - 1 -b | cos(0)|| . (13) 

The proof of the above proposition is not reported due to space limitation. 

Note that, for initial configurations satisfying (H3 we have M{xq, yo) ^ Tr ^ 
in By Proposition ^ if a reference path T has minimum radius of curvature 
Rr greater than twice the minimum turning radius R of the vehicle, then for any 
initial configuration (M{xo,yo),9o), with lateral position and orientation errors 
bounded to belong to T^~ as in , condition mu is ensured. 
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4 Verification of the Hybrid Path— Tracking Controller 



In this section the behavior of the hybrid automaton CLHA is analyzed by in- 
troducing an equivalence relation ^ in the hybrid state space O x 2?^- and by 
computing the corresponding quotient system (see jS|). 

Consider the partition U^- of the domain 2?^- g^ in (EJ given by the 

24 subsets • • • , o|, defined in Table ^ with and 

replaced by We say that {modei,{yi, 9 i)), {mode2, {y2,&2)) are 

equivalent, i.e. (modei, (yi, 6*i)) ~ (uro<ie2, (^2: ^2))) iff (yij^i) S p, for some 
p S ^{y 9)’ implies (^2)6*2) C P- We associate to the corresponding quotient 
space Q"" = {O x ■ ■ ■ ,0 x O} a nonderministic finite state machine, re- 
ferred to as FSM prp(j , with states corresponding to the equivalence classes in 
Q'^ (labeled, with a slight abuse of notation, • • • , O). The next-state func- 
tion of FSM prpfj is defined as follows: for any Qi,Q2 C Q'^, a transition from 
Qi to Q2 occurs iff there exists an arc of trajectory of the hybrid automoton 
CLHA from some {modei, (yi, 0 i)) € Qi to some {mode2, {y2, ^2)) S Q2) for some 
discrete disturbance Or and some continuous disturbance d. 

Proposition 2 . Given the hybrid system CLHA , if the discrete disturbance Ur 
takes always the value e, then, for any initial hybrid state {mode, {yo, 0 o)) € 
O X g^ as in under the action of any disturbance d bounded as in m 

with C as in JH), we have: 

— the quotient system obtained from the equivalence relation ^ is the finite 
state machine FSMpppj depicted in Figure\^ 

— an upper bound for the space travelled by the origin of the Frenet’s frame 
along the path F, when the hybrid state is in a given equivalence class is 
represented by the weight associated to exiting arc; 

— the quotient system FSMppQ remains in each equivalence class a bounded 
amount of time, except for the equivalent class O where {y, 9 ) = (0,0). 

The proof of the above proposition, which is based on reachability analysis, is 
not reported due to space limitation. 

If the reference path F has curvature always of the same sign, the convergence 
of the Dubins’ car to the path is guaranteed by: 

Corollary 1 . Lf the reference path F has curvature always of the same sign and 
amplitude lower than the hybrid feedback control and (0), ensures the 
tracking of F for any initial vehicle configuration in the domain ^{y,S) in E3). 
The origin of the Frenet’s frame covers at most a distance of 



ri + fTT+g ifce[o,^) 
(4-b77T-b^ ifC'eig^,!) 



(14) 



along the reference path F before the vehicle approaches it with correct orienta- 
tion. 
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Fig. 5. On the left: shortest paths synthesis when d = 0. On the right: quotient system 
FSM prpQ representing the behavior of the closed-loop hybrid system CLHA , when 



The proof of the above corollary is obtained by computing the longest path to 
the node O. 

By Proposition 0, if T is a straight line then the closed-loop system enforces 
sliding motions (see for a tutorial) in the space {y, 9) on the lines sr, si 
and the arcs 1^^^ until the origin is reached. If the reference path 

r is not a straight line, sliding motions are enforced only on the lines sr, si, 
on the arcs and on a piece of the arc 1*-^^ Under ideal sliding motion, 

around the origin the control uj switches at infinite frequency between ^ , 0 and 
— The mean value of such control (i.e. the equivalent control) is the signal 
kV that makes the car follows the reference path F with velocity V. In the real 
implementation smoothing techniques are applied to avoid the chattering of the 
control signal between the three values ^ , 0 and — ^ . 

The behavior of the closed-loop system CLHA under the action of the dis- 
crete disturbance Ur is characterized by the following propositions. 

Proposition 3. Given an initial condition (yo,9o) in the open neighborhood of 
the origin 

■ \y\ < 1,-arccos “ f ) < ^ < arccos Q + 0| (15) 
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Fig. 6. On the left:quotient system representing the behavior of the closed- 

loop hybrid system CLHA , when the intial state belongs to O x gy On the right: 
regions in the domain 1?^^- gj where W > Q. 



(see Figure 0), the hybrid closed-loop system CLHA keeps the continuous-time 
trajectory {y{t),9{t)) inside gy under any disturbance d{t) bounded as in 
and any sequence of events cr^.. 

Due to space limitation, the proof of the above proposition is not reported. 

Proposition 4. If the reference path F is such that changes in the curvature 
sign are at distance greater than (5 + along it, then the hybrid feedback 
control 0, with modes chosen according to m stabilizes the Dubins’ car along 
the reference path F. 

Proof. The set J\f^~ defined in (EJ is such that 

A/"(- g^ C OU U 1^^^ U U U U . 

Since, by Proposition 0 g^ is a robust invariant set for the closed-loop hy- 
brid system CLHA , then, if we restrict our attention to the domain A/j-- gy the 

transitions from lr*-^'^Ho lsr*-^^and from rl*-^'^Ho rsl^^^in the quotient system 
FSM prpQ should be removed. Furthermore, notice that, under the action of the 
discrete disturbance = switch, the reset y '■= —y and 0 := —9 introduces the 
mutual transitions iy and Hence, in the 

presence of the discrete disturbance ar and for any disturbance d as in ®, when 
the initial state belongs to O x Af^- gy the quotient system FSMppQ obtained 
from the equivalent relation ^ is as in Figure 0 

To analyse the convergence of the trajectories to O, introduce the function 

W{y,9) = ^{f + P). 



(16) 
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W{y, 9) has the property that if, at time t = t, ar = switch then W{y{i),9{t)) = 
W{y{t~),9{t~)). The derivative with respect to time of W evaluates to 



W{y,9) 



~ • o^os{0)d ~ 

ysm{9) — 0 — Ovj 

1 - yd 



V 

R' 



(17) 



where vj — and 1 in mode 2 ero, negative, and positive, respectively. The 

study of the sign of W{y, 9) is extended to the entire domain 'R(^y§y Under 
assumption (tTTl) . multiplying (HTI) by ^(1 — yd), we have 



lU > 0 O y{y, 9) = d 



—y^ sin(0) — 9cos{9) + wy9 



y sin(0) — w9 



> 0 , 



for some disturbance d bounded as in m- Hence, for any {y, 9) such that 

m{yj) = ysin(9) -079 > O, 

there exists d as in 0 such that fi{y,9) > 0 and W > 0. Otherwise, if {y,9) is 
such that r]i{y, 9) < 0, then there exists d as in 0 such that lU > 0 if and only 
if y{y,9) is positive for d = 1. That is, if 



r] 2 {y, 9) = — sm{9)y‘^ + sin(d) + w9 



9cos{9) + 079 



> 0 . 



The regions in the domain 27^ ~ where function locally increases are re- 
ported in Figure 0 Such regions are delimited by the curves rji{y,9) = 0 and 
V2{y,9) = 0. By (H3), the continuous disturbance d that maximizes W (t) is 

ifdcos(d)<0 i.e. d G (-f ,0) U (f , §7r) . . 

\o if dcos(d) > 0 i.e. d G (— |7T, — |) U (0, |) 

Consider an initial condition (yo, dp) in a neighborhood of the origin contained in 
e) Orl^^'^^. At the initial time, the hybrid model CLHA is in mode negative. 
Let us assume that Or = e, for the moment, and let us analyse the evolution 
of the hybrid model CLHA (see Figure 0. Under the action of the worst dis- 
turbance (CBI), the trajectory (y{t),9{t)) originating from {yo,9o) reaches the 
curves First W{t) decreases (in rl*-^'^^), then it increases (in rl^^’^^). Hence, 
mode switches to positive. W{t) decreases (in the first part of Ir^^'^^), and it in- 
creases again later on (in and until {y{t), 9{t)) reaches Finally, 

following a sliding motion along the curve {y{t),9{t)) reaches the origin. 

Along this trajectory W{t) assume two local maxima, which correspond to 
the intersections of 1^^^ and and two local minima: the first on the line d = 0 
when y > 0, and the second inside region Let <5 = ||(yo,do)||. Since the 

trajectory (y{t),9{t)) is continuous with respect to the initial condition (yo,dp), 
then there exist two continuous functions Cmax,Cmin : K — ^ K such that 

maxmax||(y(t),d(t))|| = Cmax(i5), minmin ||(y(t), d(t))|| = Cmin(<^)- (19) 

at at 
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Further, since the local maximum and minimum points tend to the origin as 
||(yo,^o)|| tends to zero, then lima^o Cmax(i5) = 0 and lim^^o Cmin(<5) = 0. 

Suppose now that a discrete disturbance = switch occurs at the precise 
time t at which {y{t~),9{t~)) is opposite to (yoj^o) with respect to the ori- 
gin. Then, the state is reset to {y{F),9{F)) = {—y{t~),—9(t~)) G 

■^{y 9 )^ which lies on the same line to the origin of {yo,9o). If W{ijo,9o) > 

W{y{t),9{t)) = W{—y{t~),—9{t~)) then the convergence is preserved. But, if 
W{yo,9o) < W{y{i),9{t)) = W{—y(t~),—9(t~)) then, under the action of the 
discrete disturbance ar = switch, the state (y, 9) is reset to a point farther away 
from the origin than the initial state (yg, 9 q) and convergence can be lost. 

However, if the reference path F is such that changes in the curvature sign 
are at a distance greater than {5+^)R along it, between to successive actions of 
the discrete disturbance a^, the state (y, 9) has enough time to reach the origin. 
In fact, assuming that, in the worst case, {y{F),9{t)) G fl an upper 

bound on the length the arc of F spanned by the origin of the Frenet’s frame as 
(y{t),9{t)) converges to the origin, is given by L(rl^^'^^)-|- L(rl^^'^^)-|- 
2 ^(lr(i- 2 ))_|_ that, according to the weights reported on the 

quotient system FSM prpQ depicted in Figure El evaluates to (5 + f )i?. 

To prove the robust stabilization of the car along the reference path F we have 
to show that for any e > 0, there exists <5 > 0 such that any trajectory (y(t), 9{t)) 
of the hybrid system CLHA , originating from any (yo,9o) with ||(j/o,^o)|| < <5, 
we have ||(y(t), < e- Given any e > 0, consider any initial condition {yo, 9 q) 

with 

ll(yo,0o)ll < <5 = Cix(CiUCL(e)))- (20) 

The trajectory (y(t),9(t)) evolves inside a ball of radius Cmin(Cmax(£))- If ^ 
disturbance = switch occurs at some time i, then the state is reset to 
{y{t),9{t)) = (—y{t~),—9{i~)) G §)■ In the evolution for t > t the tra- 
jectory reaches the origin before a further discrete disturbance will show up. 
Morever, since ||(y(t), 0(t)|| < Cmin(CmL(e)) then, the trajectory {y{t),9{t)) for 
t > t does not exit a ball of radius Cmax(Cmax(e)) = e- Then, the hybrid feedback 
control (0, with modes chosen according to (jOI) robustly stabilizes the car along 
the reference path F. 

5 Conclusions 

In this paper, we have used modern techniques developed for hybrid systems 
simulation and verification to solve and prove stability of a control technique 
for an interesting problem, that is route tracking by nonholonomic vehicles with 
bounds on the curvature and limited sensory information. The proposed con- 
troller is reminiscent of a synthesis proposed elsewhere for an optimal control 
problem to track straight routes, whose generalization to generic routes turned 
out to be difficult to analyze otherwise. We believe that this case study, besides 
its intrinsic interest in applications, also has a value in showing the potential of 
hybrid systems analysis techniques as applied to complex control problems. 
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